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Abstract 

Quantum corrections significantly influence the quantities observed in modern particle physics. The corre- 
sponding theoretical computations are usually quite lengthy which makes their automation mandatory. This 
review reports on the current status of automatic calculation of Feynman diagrams in particle physics. The 
most important theoretical techniques are introduced and their usefulness is demonstrated with the help of 
simple examples. A survey over frequently used programs and packages is provided, discussing their abilities 
and fields of applications. Subsequently, some powerful packages which have already been applied to impor- 
tant physical problems are described in more detail. The review closes with the discussion of a few typical 
applications for the automated computation of Feynman diagrams, addressing current physical questions like 
properties of the Z and Higgs boson, four-loop corrections to rcnormalization group functions and two-loop 
electroweak corrections. 
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1 Introduction 



1.1 Quantum field theory in particle physics 

The huge amount of data collected in the recent years mainly at the Large Electron Positron Collider (LEP) at 
CERN (Geneva) and the Stanford Linear Collider (SLC) as well as at the proton-anti-proton machine TEVA- 
TRON at FERMILAB (Chicago) has resulted in an impressive experimental precision for many parameters in 
particle physics. For example, the mass and width of the Z boson are measured with an accuracy which is of 
the same order of magnitude as for Fermi's constant PJ0]. Therefore it is essential to also improve the accuracy 
of predictions following from theoretical models. It is the purpose of this article to give an overview of some of 
the technical tools to perform the corresponding calculations. 



ai piesein, int; must lavuuieu iiieuieiiccii muuei 111 panicie piiysics is me su-cciiieu oiaiiuaru iviuuei. 11 nicui- 
porates all known quarks and leptons and describes their interactions by the gauge group SU(3) xSU(2) xU(l). 
Here, the SU(3) generates the strong interactions among the quarks by the exchange of eight massless gluons, 
while the SU(2)xU(l), usually called the electroweak sector, comprises electromagnetic and weak interactions, 
carried by the massless photon and the massive charged and neutral Z boson, respectively. 

Up to now all theoretical predictions deduced from the Standard Model have found perfect confirmation in 
experiment. For example, the existence of the top quark was predicted after the discovery of the T-resonances 
in 1977 H merely out of the need for an SU(2) doublet partner for the bottom quark. At least equally 
impressive, however, was the accurate prediction of its mass, Mt, by comparing precision measurements of 
certain observables at LEP and SLC to their theoretical values. The top quark mass dependence in these 
observables only appears as a higher order effect through the so-called radiative corrections, typical for quantum 
field theories. So when the first evidence of the discovery of the top quark was announced by the CDF and DO 
collaborations the agreement of directly and indirectly measured values for the top mass served as a very 
convincing argument. Today statistics have improved considerably and the original signal has been impressively 
confirmed ||. 

A similar situation exists with the Higgs boson, at the present time being the only particle of the Standard 
Model not yet discovered. Its existence, however, is by no means guaranteed as was the top quark after b 
discovery. Furthermore, predictions for the Higgs mass Mh suffer from much larger uncertainties because 
the dependence of the relevant radiative corrections on Mjj is only logarithmic || while it is quadratic on Mt. 
Nevertheless, it is possible to set the limits to Mh = 76^^ GeV with an upper bound of Mh < 262 GeV at 95% 
confidence level [jTj. The current lower limit for Mh from the direct search at LEP reads M#;>89.8 GeV ||. 

From the theoretical point of view the Standard Model does not look like a "final" theory. First of all, 
it contains a lot of free parameters and any attempt to predict their values needs a generalization of the 
model. Furthermore, there is no way to find unification of the strong, weak and electromagnetic coupling 
constant within the Standard Model. A very promising and theoretically attractive extension is the Minimal 
Supersymmetric Standard Model (MSSM) which not only allows for gauge unification but even solves the 
problem of quadratic divergences and many others. Phenomenological calculations, however, are even more 
difficult here since the particle content is approximately doubled. Moreover, some important technical tools 
are no longer applicable. 

Let us now turn to the question how to extract physical information like cross sections and decay rates from 
the Standard Model, leading to predictions like the values for Mt and Mh mentioned above. At the moment 
the vast majority of the calculations are based on perturbation theory, represented by the extremely depictive 
Feynman diagrams. Together with the Feynman rules, which constitute the translation rules from the graphical 
to a mathematical notation, they are in complete one-to-one correspondence to the individual terms of the 
perturbative series. As long as one is interested only in perturbative results, one may formulate the theory in 
terms of the Feynman rules alone, i.e., by providing the list of possible vertices, describing particle interactions, 
and propagators. Any physical process may then be evaluated by connecting the vertices and propagators in all 
possible ways to give the desired initial and final states. The specific order in perturbation theory is reflected 
in the number of vertices of the corresponding Feynman diagram. Given a fixed initial state and increasing the 
number of vertices leads to the appearance either of closed loops or of additional particles in the final state. 
In the former case, this implies an integration over the corresponding loop momentum, while the latter case 
results in a more complicated phase space integral when calculating the rate. Usually, diagrams with more 
than one closed loop are called multi-loop, those with more than two particles ("legs") in the final state are 
called multi-leg diagrams. The technology, in particular the automation of the computations in the multi-loop 
sector is much further developed than for multi-leg diagrams. Furthermore, to a certain extent it is possible to 
relate multi-leg diagrams to multi-loop ones via the optical theorem. Thus our main focus in this article lies 
in the calculation of loop diagrams. 

1.2 Motivation for automatic computation of Feynman diagrams 

The first realization of the idea to pass purely algebraic or, more precisely, symbolic operations to a computer 
indeed was driven by particle physics, when in 1967 M. Veltman developed the program SCHDDNSCHIP Q, 
mainly to control the evaluation of fermion traces, i.e. traces of gamma matrices. Today, SCHOONSCHIP to a large 
extent is superceded by other systems which have, however, taken over many of its basic ideas. Some important 
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others. The by now highly developed algebraic system Mathematica |]l2[ actually is a derivative of the program 
SMP, created in 1980 by the particle physicist S. Wolfram. While some of these programs claim to be fairly 



general, the most direct descendant of SCHDDNSCHIP is FORM [13], a system mainly tailored to high energy 



physics. In Section a few representatives of the most important algebraic systems will be introduced. The 
main concern of this review are not these algebraic programs, but rather their applications, i.e., packages based 
on these systems. It should be clear, however, that they are very important prerequisites for the automation 
of Feynman diagram calculations. 

The higher the order of perturbation theory under consideration, the larger is the number of contributing 
Feynman diagrams and the more complicated is their individual evaluation. This makes it desirable and even 
unavoidable to develop efficient algorithms which can be implemented on a computer in a simple way. Therefore, 
large efforts in the field of Feynman diagram evaluation are devoted to find algorithms for the specific types 
of operations arising in a typical calculation. One may characterize these operations by analyzing the steps 
necessary for a perturbative field theoretic calculation: given a set of Feynman rules, first the contributing 
Feynman diagrams for the process under consideration have to be generated. This, being mainly a question 
of combinatorics, is clearly a well suited problem for computerization, and meanwhile very effective algorithms 
have been implemented with emphasis on slightly diverse purposes. The next step is the evaluation of the 
corresponding mathematical expressions for the diagrams which can be divided into an algebraic and an analytic 
part. The algebraic part consists of operations like contracting Lorentz indices and calculating fermion traces. 
Nowadays many algebraic systems provide functions optimized for this kind of manipulations. 

The analytic part of a Feynman diagram calculation is concerned with the Feynman integrals, and usually this 
is the point where computers fail to be applicable without human intervention. However, Feynman integrals 
establish a very special class of integrals, and for huge subclasses algorithms exist that algebraically reduce each 
of them to basic sets of meanwhile tabulated expressions. Two of the most important techniques of this kind 
are the so-called tensor reduction- [14] and the integration-by-parts algorithm [15|. The first one reduces any 



Lorentz structure in the integrand to invariants and has been worked out explicitly at the one-loop level for 
the general case and up to two loops for propagator diagrams. Using the integration-by-parts algorithm, on 
the other hand, recurrence relations can be derived which express diagrams with non-trivial topologies through 
simpler ones at the cost of increasing their number and the degree of the denominators. The practical realization 
of these algorithms would be unfeasible in realistic calculations without the use of powerful computer systems. 

From the considerations of the previous paragraphs one can see that as long as one is interested in the exact 
evaluation of Feynman diagrams, the number of problems that can completely be passed to computers is rather 
limited which is mainly due to the relatively small number of exactly solvable Feynman integrals. There are, 
however, many processes where some kind of approximate result may be equally helpful. This becomes clear 
by recalling that, on the one hand, working at fixed order in perturbation theory is an approximation anyway, 
and, on the other hand, experimental results are not free of uncertainties either. Therefore, it is enormously 
important to have approximation procedures for Feynman diagrams. For instance, the numerical evaluation of 
complicated Feynman integrals, which is a typical task for a computer, is a possibility to obtain results with finite 
accuracy. This approach is currently favoured for the automation of tree- and one-loop calculations involving 
lots of different mass parameters like in the Standard Model or even the MSSM. For recent developments in 



the numerical evaluation of Feynman diagrams beyond the one- loop level we refer to [16] 



A different kind of approximation may be obtained by expanding the integrals with respect to ratios of different 
scale parameters. Here the so-called asymptotic expansions of Feynman diagrams are becoming more and more 
popular, since they allow any multi-scale diagram to be reduced to single-scale ones, at least in principle. While 
the evaluation of the former may be a hard job even at one-loop level, the latter ones can be computed up 
to three loops using the integration-by-parts algorithm mentioned before. The required efforts for a manual 
application of asymptotic expansions increase steeply with the number of loops. Meanwhile, however, the 
most important variants have been implemented and have also been applied to the calculation of three-loop 
contributions to several important physical observables. 

Performing huge calculations with computers inevitably leads to seemingly unrelated difficulties of purely 
organisational nature: Computers nowadays are not yet as stable as one would like them to be. In this way 
they force the user to cover on possible breakdowns, keeping the loss minimal. In the case where the long 
run-time is mainly due to a huge number of diagrams, each of them taking only a few minutes to be evaluated, 
the solution to this problem is rather straightforward. Saving each diagram on disk after computation naturally 
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the like). For more complicated diagrams, however, it is helpful if the algebra program itself is structured in a 
way that intermediate results will not be lost after system crashes. 

All the algorithms mentioned above allow the calculation of Feynman diagrams to be automated to more or 
less high degree, where one will certainly choose different combinations of them for different purposes. It is 
quite difficult to give an exhaustive definition of automation. For instance, often the computer system either 
applies only for one special problem (or a very restricted class of problems), nevertheless saving the user a 
great deal of effort and preventing otherwise unavoidable errors. In other cases the system is rather flexible 
but requires a certain amount of human intervention for each specific problem. The final goal of automation 
where a typical screen snapshot would look like 

model: SM 
loops: 3 

initial-state: ep,em 
final-state: t,tb 
parameters: default 
observable: sigma_tot 
energy: 500 GeV 

sigma = 0.8 pb 

is at the moment still out of range. Nevertheless, many of the recently performed calculations in high energy 
physics would never have been possible without a high degree of automation. Let us, at the end of this 
introduction, list a few of these applications that will be addressed in more detail in Section ^[ 

1.3 Examples for physical applications 

Most of the multi-loop and multi-leg computations have been performed for the QCD sector of the Standard 
Model. On the one hand they are numerically most important simply because the coupling constant is relatively 
large. On the other hand they are simpler to evaluate due to the fact that very often only one dimensionful 
scale appears like the external momentum or a large internal mass. 

The "classical" example which is currently known to order are QCD corrections to the hadronic R ratio, 
R(s) = a(e + e~ — > hadrons)/cr(e + e - — > fi + fi~). Already in 1979 0(a 2 s ) corrections have been computed [17 



in the massless limit. A large part of the calculations was carried out without using computers and lots of 
efforts had to be spent on nowadays trivial things like fermion traces. It took more than ten years until the 
next order was completed [18|, Is|| . Even then the manipulation of the diagrams was largely done "by hand". 



As a consequence, important checks like gauge parameter independence could not be performed. Only recently 
the O(af) result was checked by a completely independent calculation which used the powerful methods of 
automatic generation and computation of diagrams (T^j . 

The full mass dependence of the 0{a 2 ) corrections to R(s) has become available |2(| only very recently, when 
high moments of the polarization function were evaluated. In this calculation the problem was not in the 
number of diagrams which is less than 20, but in the huge expressions occurring in intermediate steps. (For 
some diagrams up to four giga- bytes of disk space were required.) For such computations it is indispensible 
that enormously powerful algebra systems are available and that their use is completely automated to avoid 
any possible human error in intermediate steps. 

The number of diagrams to be considered for the 0(af) corrections to R(s) at m = is of the order of one 
hundred and thus it is still possible to organize the calculation by hand. The situation is different for the decay 
of the Higgs boson into gluons, for example. Among other reasons this is a very important process because 
the inverse reaction of gluon fusion is the dominant production mechanism of Higgs bosons at the future Large 
Hadron Collider (LHC). Since the leading order QCD corrections to the loop-induced process H — > gg amount 



to roughly 68% |pl, 22|, it was necessary to evaluate the next-to- leading order contribution. However, the 
number of three-loop diagrams to be computed in this case is of the order one thousand. This in combination 
with the large expressions in the intermediate calculational steps is the reason why one has to rely on a high 



degree of automation as will be described in Section 4.4 
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governing the running of the strong coupling constant and the quark masses, respectively, were computed at 
four-loop level. Roughly 50 000 diagrams for (3 Q] and 2 000 for j m |?4|] contributed, and it is quite clear 
that such calculations would never have been possible without the intensive use of extremely efficient computer 
algebra systems in combination with administrative software concerned with book-keeping. 

So far, none of the examples above really had to rely on one of the approximation procedures like numerical 
integration or asymptotic expansion mentioned above. This becomes relevant for the decay of the Z boson 
into quarks at 0(aa s ), for example. Here, the bottom quark channel is of special interest since the top quark 
appears as a virtual particle. With the present technology, a calculation of the full M^-dependence is out of 
reach because up to four different mass scales are involved. However, asymptotic expansions provide a very 
promising method to get a result which is almost equivalent to the full answer. In the approach which was used 



to tackle the 0(aa s ) terms to this process [25], the 69 contributing diagrams are split into 234 subdiagrams. 
The manual application of the method of asymptotic expansions was therefore impossible and the existence of 
program packages that apply this procedure automatically was important. 

Meanwhile, experimental accuracy allows the electroweak sector of the Standard Model to be tested even at 
the two-loop level, and several groups have tackled the corresponding calculations. As was already outlined, 
the strategy here is to reduce all integrals to scalar ones. The latter depend on a lot of different mass scales, 



such that their evaluation must be done numerically. In Section 4.6 we will describe the computation pq] of 



two- loop corrections induced by the Higgs boson to Ar which enters the relation between Gf, Mz, M\y and a. 

It is not only since the increase of the center-of-mass energy at LEP above the W boson production threshold 
that processes involving four (or more) particles in the final state became very topical. In such reactions even 
the contribution from lowest order perturbation theory may pose a serious problem as quite a lot of diagrams 



are involved and a highly non-trivial phase space integration has to be performed. In Section 4.5 the scattering 



of vector bosons is discussed in the background of a potentially large coupling among the bosons at high 



energies [27|. 



Other applications that will not be discussed in more detail but shall further substantiate the success of 



computer algebra in Feynman diagram calculations are, for example, QCD corrections to Ap [28, 2{| and Ar [3C] 



moments of QCD structure functions [31-33], top mass effects in the decay of intermediate [p4l and heavy 



Higgs bosons as well as in e + e~ collisions [36|. Involved one- loop calculations concerning radiative corrections 
to the gauge boson scattering were considered in fl3"7| . The muon anomalous magnetic moment [|38| and the 
neutron anomalous electric moment |3{| are important examples of two-loop calculations in the electroweak 
sector. In the main difficulty was the evaluation of two-loop vertex diagrams in order to obtain corrections 



to the decay of a heavy Higgs boson. Investigations of b decay [41, 42] and threshold production of heavy 
quarks p3|-446| are examples for the developments of new algorithms and their implementation. A large list 
of results obtained with the help of programs that automatically compute tree-level processes involving many 
particles can be found in |l7]] . 

Many more important calculations for physical observables in the multi-loop and multi-leg sector have been 
performed — mainly within the last few years, but let us close the list at this point. 

The outline of this review is as follows: Section ^ introduces the most important technical tools that allow the 
automation of Feynman diagram calculations. This includes on the one hand some of the required algorithms, 
and on the other hand the most frequently used algebraic programming languages to implement these algo- 
rithms. Concrete examples for such implementations will be described in Section |3|. This splits into a surveying 
part whose main purpose is to collect a list of the most important programs, and a more specific part concerned 
with the actual realization in the light of a representative set of selected packages. Finally, Section |] gives a 
list of physical applications that have been considered in the literature and that would not have been feasible 
without the help of a certain degree of automation. Some of these examples actually were the driving force for 
some of the afore-mentioned packages to be developed. Others are real applications in the sense that existing 
programs could successfully be used to tackle as yet unsolved problems even by people that were not involved 
in the development of the programs. 



2 Theoretical tools 
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The main field for multi-loop calculations certainly is QCD in the perturbative regime. The reason is, on the 
one hand, that in general the complexity of a given Feynman diagram increases rapidly with the number of 
dimensional parameters involved. In the Standard Model, for example, with its relatively large amount of 
different particles and scales, with the currently available technology it is hardly possible to compute processes 
for electroweak phenomena beyond two-loop level. In QCD, by contrast, the gauge bosons are massless, and 
the quark masses are such that for many physical reactions one may either neglect them completely or consider 
only one quark flavour as massive and all others as massless particles. On the other hand, in QED, the fine 
structure constant is very small, rendering the higher order corrections negligible in general. The coupling 
constant in QCD is about ten times larger, but still small enough to play the role of an expansion parameter. 

2.1.1 Dimensional Regularization and Minimal Subtraction 

The momentum integrations of loop diagrams are in general not convergent in four-dimensional space-time. 
This requires the introduction of the concept of renormalization, rendering physical quantities finite and at- 
tributing a natural interpretation to the parameters and fields of the underlying theory. In order to isolate 
the divergent pieces one introduces a so-called regularization scheme, the most popular one being dimensional 
regularization E3| at that time. The strategy is to replace the four-dimensional loop integrals appearing in 
Feynman diagrams by "D-dimensional" ones, obeying the basic relations of standard convergent integrals like 
linearity etc. D is considered to be a complex parameter. The limit of integer D, if it exists, is required to 
reproduce the value for the standard integral. All manipulations are those of convergent integrals, and only 
after integration one takes the limit D — ► 4. The divergences are then reflected as poles in 1/(D — 4). 

The concept of renormalization consists in compensating these divergences by adding suitable divergent terms, 
so-called counter-terms, to the Lagrange density. Since the only requirement on the counter-terms is to cancel 
the divergences, there obviously is a freedom in choosing their finite parts. A particular choice is called a 
renormalization scheme. 

A convenient renormalization scheme in combination with dimensional regularization is the so-called minimal 
subtraction or MS scheme JKJ. It prescribes to precisely subtract the poles in 

e = (4-D)/2. (1) 

Even more popular is the MS scheme (modified MS scheme) |>(| which is based on the observation that the 
poles in 1/e appear only in the combination 

1 

A = --7E + ln4vr (2) 

and that therefore it is convenient to subtract A instead of 1/e. In the following we will use the MS scheme 
throughout unless stated otherwise. If divergent parts will be quoted explicitly, the 7e and In Air terms will be 
omitted. 

Dimensional regularization, like any other regularization scheme, introduces an arbitrary mass parameter, 
usually denoted by [i. This becomes clear by considering, for example, the one-loop integral 

1 f d D p 1 _ {-l/q 2 ) a+b - D l 2 T{a + b- D/2)Y{D/2 - a)T{D/2 - b) 

~ij (27r)^(_( g -p)2)a(_p2 ) 6 (tojm T(D -a- b)T(a)T(b) ' ^ 

where T(x) is Euler's gamma function. Note that a and b are not necessarily integer. The Laurent-series of the 
r.h.s. with respect to e leads to logarithms with dimensional arguments. Therefore, one multiplies any quantity 
by an appropriate D-dependent power of fi such that the mass dimension of the whole object is independent 
of D. It is clear that only logarithmic /^-dependence may arise by this procedure. In terms of the Lagrangian, 
the artificial mass scale manifests itself in a coupling constant with D-dependent mass dimension. 

The combination of dimensional regularization and minimal subtraction has many computational consequences. 
One can show, for example, that within the framework of dimensional regularization massless tadpoles, i.e. in- 
tegrals that do not carry any dimensional parameter except the integration momenta, may be set to zero 
consistently. On the other hand, minimal subtraction guarantees that any renormalization constant is a series 



in nit; coupling consiam aione, wiiiioui explicit uepenuence on any uimensioiiai quamiiy nite masses ur mo- 
menta. Many calculations can be considerably simplified by exploiting one of these properties. Indeed, one 
of the most powerful tools for the computation of multi-loop diagrams, the so-called algorithm of integration- 
by-parts, strongly resides on the properties of dimensional regularization. It will be described in more detail 
below. 

There are also certain drawbacks of dimensional regularization as well as the MS scheme. One of them, of 
course, is the lack of any reference to physical intuition, as one has it for regularization schemes like introducing 
cut-offs or a discrete space-time lattice, and for renormalization schemes like the on-shell scheme. Besides this, 
dimensional regularization generally causes problems whenever explicit reference to four dimensions is made. 
For example, the anti-commuting definition of 75 leads to inconsistencies when working in D dimensions. A 



consistent definition of 75 was given in [Eq| and formalized in [51|. It defines 75 to anti-commute with 70, . . . , 73 
and to commute with all the remaining 7-matrices. This definition obviously breaks Lorentz covariance and 
requires the introduction of so-called evanescent operators which makes practical calculations quite tedious. 
Also in the Supersymmetric Standard Model, defined in four-dimensional space time, one runs into problems, 
but we will not dwell on them here since this will not concern what follows. 



2.1.2 Integration-by-parts 

Let us now describe one of the benefits of dimensional regularization in more detail, namely the integration- 
by-parts algorithm. It uses the fact that the Z)-dimensional integral over a total derivative is equal to zero: 

/d^A /(Pj ...) = . (4) 

By explicitly performing the differentiations one obtains recurrence relations connecting a complicated Feynman 
integral to several simpler ones. The proper combination of different recurrence relations allows any Feynman 
integral (at least single-scale ones) to be reduced to a small set of so-called master integrals. The latter ones 
have to be evaluated only once and for all, either analytically or numerically. 

The integration-by-parts algorithm was initially introduced for massless two-point functions up to three loops 



15 1, where two non-trivial master integrals were to be evaluated. Further on, the technique was applied to 
those massive tadpole integrals contributing to the three-loop QCD corrections of the photon polarization 
function in the limit q 2 <C m? [52|, where q is the external momentum of the correlator and m is the mass of 



the heavy quark. One non-trivial master integral results in this case. The procedure was extended to apply 



to the three- loop QCD corrections for the p parameter [28, 29| where two more master integrals had to be 
evaluated. 



The recurrence relations for all possible three-loop tadpole integrals with a single mass were derived in |5J 
and the (most complicated) master integrals were calculated in 54]. At four- loop level the integration- by-parts 



technique was applied to completely massive tadpole diagrams, only aiming at their divergent parts, however. 
This restricted problem leads to two four-loop master integrals p3[| . 

The integration- by-parts technique equally well applies to on-shell integrals, the complexity at n-loop level 
being comparable to the tadpole case at n + 1 loops, however. The recurrence relations at two loops were 



derived some time ago [55| and were applied to the fermion propagator in order to determine the relation 
between the on-shell and MS mass in QCD ]5q| , and the wave function renormalization constant |56[ . The 
three-loop on-shell integrals contributing to the anomalous magnetic moment of the electron could be reduced 
to 18 master integrals with the help of the integration-by-parts algorithm, and thus an analytic evaluation of 
this quantity could be performed ^7|. Very recently 0(a 2 ) corrections of the \i decay were calculated f58j| , and 
the integration-by-parts method was used to determine the pole part of the corresponding four-loop on-shell 
integrals. 

To demonstrate the power of the integration-by-parts method let us consider the scalar two-loop diagram of 
Fig. ||. The corresponding Feynman integral shall be denoted by 

f d D p d D k 1 
/(ni, . . . , n 5 ) = J (27r)D (27r)D ^ + TO 2)m . . . (pg + m |)n 5 ' ( 5 ) 



where Pi, ■ ■ ■ ,Ps are combinations of the loop momenta p, k and the external momentum q (we work in Euclidean 
space here), m, . . . ,n$ are called the indices of the integral. Consider the subloop defined by lines 2, 3 and 5, 




Figure 1: Two-loop master diagram. The arrows denote the direction of momentum 
flow. 



and take its loop momentum to be p = p$. If we then apply the operator (d/dp§) -p?, to the integrand of I, we 
obtain a relation of the form (0), where 



/(P5,...) 



(pi + m\Y^p\ + m\Y^p\ + m\) 



«3 



(6) 



Performing the differentiation and using momentum conservation at each vertex one derives the following 
equation: 



n 3 



3+ ({ 



9 9 9 

4 + m 4 — m 5 — m 3 



n 2 2 + (5 — 1 + ra\ — m\ — m\ 

2 



+D - 2n 5 - n 3 - n 2 + 2n 5 m 5 5 + J(rai, . . . , n 5 ) 



0. 



(7) 



where the operators l 1 * 1 , 2 , . . . are used in order to raise and lower the indices: I ± /(. . . , rii, . . .) = /(..., rii ± 
1, . . .). In Eq. (|7|), generally referred to as the triangle rule, it is understood that the operators to the left of 
I(n\, . . . ,715) are applied before integration. If the condition 7775 = 0,7713 = 7774 and mi = 777,2 holds, increasing 
one index always means to reduce another one. Therefore this recurrence relation may be used to shift the 
indices ri\, 774 or 775 to zero which leads to much simpler integrals. 

The triangle rule constitutes an important building block for the general recurrence relations. The strategy 
is to combine several independent equations of the kind (|7|) in order to arrive at relations connecting one 
complicated integral to a set of simpler ones. For example, while the direct evaluation of even the completely 
massless case for the diagram in Fig. [I] is non-trivial, application of the triangle rule (0) leads to 



I(ni, 



1 



,n 5 ) 



D - 2n 5 - 77 2 - 77 3 



n 2 2- 



I") + 77 3 3+ (5- 



J(m,. 



.n 5 ) 



(8) 



Repeated application of this equation reduces one of the indices 771, 774 or 775 to zero. For example, for the 
simplest case {n\ = n 2 = . . . = 775 = 1) one obtains the equation pictured in Fig. |2j: The non-trivial diagram 
on the l.h.s. is expressed as a sum of two quite simple integrals which can be solved by applying the one-loop 
formula of Eq. (^) twice. This example also shows a possible trap of the integration-by-parts technique. In 
general its application introduces artificial 1/e poles which cancel only after combining all terms. They require 
the expansion of the individual terms up to sufficiently high powers in e in order to obtain, for example, the 
finite part of the original diagram. This point must carefully be respected in computer realizations of the 
integration-by-parts algorithm: One must not cut the series at too low powers because then the result goes 
wrong; keeping too many terms, on the other hand, may intolerably slow down the performance. 

In our example, the l.h.s. in Fig. [2] is finite, each term on the r.h.s., however, develops 1/e 2 poles. The first 
three orders in the expansion for e — > cancel, and the 0{e) term of the square bracket, together with the 1/e 
in front of it, leads to the well-known result: 1(1, 1, 1, 1, 1) = 6C(3)/g 2 , where q is the external momentum. 




<x> 



Figure 2: Symbolic equation resulting from Eq. (||) applied to the diagram 
1(1, 1, 1, 1, 1). The dot indicates that the respective denominator appears twice. 



in general, me successive appiicauun ui recurreuw; reiauuns geneiaues a imge numuer ui lerms uui ui a single 
diagram. Therefore, a calculation carried out by hand becomes very tedious and the use of computer algebra 
is essential. 

At the end of this section let us describe an alternative approach which tries to avoid the explicit use of 
recurrence relations and thus the large number of terms in intermediate steps of the calculation. The crucial 
observation is that an arbitrary integral is expressible as linear combination of the master integrals where the 
coefficients simply depend on the dimension D and the indices of the original integral. Therefore, an attempt 
to explicitly solve the system of recurrence relations in terms of integral representations was made in |5^,|6(| . It 
was even possible to derive additional recurrence relations over the space-time dimension D in this approach. 
It was successfully applied to the class of three-loop tadpoles pictured in Fig. |3[ where a significant reduction 
of CPU time could be achieved. Further developments in this direction look quite promising. 




Figure 3: Three- loop topology for which the system of recurrence relations was 
explicitely solved. Solid lines carry a common mass M, dashed lines are massless. 



2.1.3 Tensor decomposition and tensor reduction 

The calculation of Feynman diagrams for realistic field theories inevitably leads to tensor integrals, i.e. Feynman 
integrals carrying loop momenta with free Lorentz indices in the numerator. Since allowing for tensor structure 
largely increases the number of possible integrals, it is important to have algorithms that reduce them to a 
set of basis integrals. The algorithm described in the present section is referred to as the Passarino-Veltman 
method |14}|. At one loop-level, it reduces any tensor integral to integrals with unity in the numerator. At two- 



loop level this is no longer true, except for two-loop propagator-type integrals [61|. We will briefly introduce 
the corresponding technique at the end of this section. 

Consider an arbitrary 1-loop integral carrying tensor structure in the integrand, 



I (27r) D DqDi ... Dn _ 1 -\ DqDi • • • Djv-1 / D ' M 



with 

D = k 2 - ml D i = (k+ Pi ) 2 - m 2 (10) 

(the ie in the denominator is suppressed here). The statement is (l4j that it can be expressed in terms of the 
set of scalar integrals defined as 

I N = U n 1 n ) (U) 



D D 1 ■■■d n . 



1/ D 



The strategy is more or less straightforward and relies on the decomposition of the integral under consideration 
into covariants, built out of g^ u and the external momenta. The general solution to this problem in terms of a 
recursive algorithm can be found in |^]. It reduces the rank of the tensor in the numerator by one in each step. 
Here we only want to give an example to get an idea of how the algorithm works. Consider the three-point 
second-rank tensor integral 

1 = (i££n:) . ( 12 ) 



D Q DiD 2 / D ' 



with its tensor decomposition 



tfcv = COOSV + CllPljuPlf + C 2 2P2fiP2u + Ci 2 Plfj,P2u + C 2 lP2^Plu ■ (13) 



j_ne Ojj are uneii caneu russurniu- vtiLiiiun cuemcienis aim are iiameu uy /i, cuirespunumg iu uie 

value of JV = 1, 2, 3, ... in Eq. (|). Denning 



Rt 



k»( Pi -k)\ 



and 



we may write the tensor decomposition of those quantities as 

R? = ni pt + r i2 pti, 



R 



k 2 



00 



D D 1 D 2 



(14) 



Contracting (|13|) by 
Pi and ^2 : 



00 = ^00 • (15) 
i = 1,2) one obtains a set of four equations by separately comparing coefficients of 

2 



rn = Coo + Cnpf + Ci2Pi-P2, 

T12 = C 2 2Pl'P2 + C21P1, 

+ C 2 1P1-P2, 



T22 
T21 



Coo + C22 p\ 



(16) 



Cupi-P2 + C12P2. 



Contraction of ( |13| ) with g^, on the other hand, yields 

^00 



DC 00 + Cnpl + C22P2 + C12P1 -P2 + C21P1 -P2 ■ 



(17) 



Together, these are five equations for the five unknowns Cy (actually there are only four unknowns, since 
C12 = C21 because of the symmetry \i <-> v in 1^; this redundancy may serve as a useful check in the end). 
Combining (|17]) with the sum of the first and third equation in ([H]), one immediately gets the solution for Coo: 
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(roo - rn - r 22 ) 



(18) 



The remaining coefficients may be determined by rewriting (^) as two sets of systems of linear equations: 



ru - Coo 
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C11 
C12 



where 
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Pi -P2 



T\2 
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Pi -P2 
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C21 
C22 



(19) 



(20) 



So, if X is invertible, Eq. ( |i~8| ) and the inverse of Eqs. ( |l~9| ) determine the Cij in terms of the ry. 
In turn, by rewriting 

k-pi = -[Di - D - fi] , with /j = p, 2 - mf + ml , 
and inserting this into the first equation of (|i~4]), one arrives at 
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(22) 



which allows to compute the through first-rank tensor integrals. Thus, the first run-through of recurrence 
is done. The second one, in turn, is only concerned with at most first-rank tensors. In that way, any tensor 
integral can be reduced to the basic set of scalar integrals defined in (11), provided that X is invertible. Several 
improvements to the original algorithm concerned with the problem of vanishing det X and also with the 
numerical evaluation of the scalar integrals have been worked out An overview and a complete list of 

references can be found in 1651 . 



At two-loop level this strategy no longer reduces numerators of the integrals to unity. The reason is that 
in general one encounters "irreducible numerators", i.e. scalar products that are not expressible in terms of 
the denominator via relations like (^l|). But in the case of propagator-type integrals, i.e. those with only one 
external momentum p, the full reduction may still be achieved by applying the mechanism above first to a 
subloop |pl). For example, consider the diagram in Fig. 0, and a corresponding first-rank tensor-integral: 




Figure 4: Two-loop propagator-type diagram. The momenta ki are linear combina- 
tions of the loop momenta I, k and the external momentum p. 



S = [ d ° k d ° l (23) 

The notation here is a bit different from the one in (|l0|): 

D ki = kf-m 2 , (24) 

with 

&)2 = k , &3 = k + 1 , &4 = / + p , k§ = 1 . (25) 

Direct application of the Passarino-Veltman algorithm through the ansatz Sn = p^S(p 2 ) and contracting this 
equation by p^ leads to the irreducible numerator p-k, so that no simplification results in this way. 

The idea is instead to perform the tensor reduction of the subintegral over k first, considering I as its external 
momentum: 

c = f d n l _l_f d D k k^ = f d D l _^.l 4 f d D k . 



" J DiDp+iJ D k D k+ i J DtDp+i I 2 J D k D k+l 



l 2 DiD p+ J D k D k+ i p 2 

p v ( aDi. f aDi [gfc+j ~ Dk-l 2 ~m 2 + ml] [D p+l -D l -p 2 -m 2 5 + mj] 
-^J dk J dl l 2 D k D k+l D lDp+l • (26) 

After canceling common factors, the numerator has no dependence on loop momenta any more which is what 
we were aiming for. The application to arbitrary propagator-type tensor integrals can be found in 



2.1.4 Tensor reduction by shifting the space-time dimension 



As already noted in Section 2.1.2, the Passarino-Veltman method reduces the tensor structure in the numerator 
to unity only at one-loop level, the reason being "irreducible numerators" appearing at higher loop order. An 
alternative approach that circumvents this problem has been worked out at one-loop level in |3(| and generalized 
to an arbitrary number of loops in [p^,^8|. However, at two-loop level it has only been applied to propagator 
type integrals up to now (see, e.g., |5J|). The basic idea is to express tensor integrals in D dimensions through 
scalar ones with a shifted value of D. Again we do not want to present the algorithm in its full generality here, 
but try to shed some light on the main ideas by giving a concrete example. 



Consider again the integral (|l2|). By introducing an auxiliary vector a^, it may be written as 
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Using the Schwinger-parameterization for propagators 
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(29) 



where da = daodaida 2 , and with the help of 

d D k 



(2k) 



D 



exp 



i(Ak 2 + 2p-k) 



e 1 a 



' (4vri A) D / 2 ' 



(30) 



one obtains 
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D DiD 2 ~ l D (4iri) D / 2 J (a + a x + a 2 ) D / 2 
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Inserting (|3l|) into (pTj) and explicitly performing the differentiations yields 



ai(a + a 2 )p\ 



(31) 
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(32) 



The technique to derive Eq. (^) and its generalization to arbitrary multi-loop diagrams is known since long 
(see, e.g., []5l|). However, instead of differentiating with respect to a we may equally well apply the following 
operator on the scalar integral ^ £)oD 1 D 2 ) D ' m wa ^ Siting rid of the auxiliary vector a: 

T d \ n d \ 

exp [-i(a\pi ■ a + a 2 p 2 ■ a)p\ 



a=0, oij=i— Q-tt, p=4-7rid+ 



d d \ ( d d , 
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(33) 



where the operator d + increases the space-time dimension by 2, i.e. d + (• • -) D = (■ ■ -)d +2 - Finally, we have 
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(34) 



The algorithm formally applies to an arbitrary number of loops and external legs. But this means only that 
any tensor integral can be reduced to scalar integrals with a shifted number of space-time dimension, the latter 
ones remaining still to be evaluated. A strategy to cope with these diagrams is to use generalized recurrence 
relations fl68fl . At the one-loop level they have been worked out for arbitrary n-point functions in [p^| . At the 
two-loop level, however, so far they are only published for propagator-type diagrams |68|| . Since we feel that 
these generalized recurrence relations are beyond the scope of this review, let us refer the interested reader to 

The algorithm described above was used in |39| for the computation of the two-loop 
QCD corrections to the fermion propagator. 



the literature 
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2.2 Asymptotic expansion of Feynman diagrams 



a.a.x. dentil aiiLies 



As was outlined in Section the integration-by-parts algorithm was successfully applied to single-scale 
integrals, i.e. massive tadpole, massless propagator-type, or on-shell integrals. For an arbitrary multi-scale 
diagram it is in general rather difficult to solve recurrence relations. 

However, if the scales involved follow a certain hierarchy, a factorization is possible. For example, consider the 
operator product expansion for the correlator of currents j(x) = ^(x)Txp(x) in the limit Q 2 —* oo, where T is 
some Dirac matrix and ip a quark field with mass m: 

i f dxe^ x Tj(x)j(0) Q2 ~°° £ C n O n , (35) 



where T denotes the time ordered product. It was realized in [71-74 that if one adopts the minimal subtraction 
regularization scheme and abandons normal ordering of the operators, then Q appears only in the coefficient 
functions C n while the mass m is completely absorbed into the operators O n and appears only in the matrix 
elements. 

The attempts to find an algorithm that produces this factorization for arbitrary Feynman integrals finally 
resulted in the prescriptions for the asymptotic expansion of Feynman diagrams ]75|-[79[| (see also |)80|| ). These 
prescriptions provide well defined recipes that are completely decoupled from any field theoretic derivation and 
even lack a rigorous proof, as in the cases (iii) and (iv) below. However, their success in practical applications 
justifies them a posteriori. 

At the moment the following procedures are used in the calculations: 



(i) Large- Momentum Procedure: Q 3> q,m 

(ii) Hard-Mass Procedure: M 3> q,m 

(iii) Threshold Expansion 

(iv) Expansion with the external momenta on the mass shell. 



The first two of them will be considered more closely in the next subsection. The presentation will be rather 
informal, explaining the procedures in a ready-to-use form. Currently, the technical apparatus for the latter 
two cases is much less developed. For this reason they will be touched upon only briefly here. 

One may treat large-momentum and hard-mass procedure on the same footing. Thus, in what follows we only 
present the general formulae in the case of large external momenta — the transition to the hard-mass procedure 
is straightforward. The prescription for the large-momentum procedure is summarized by the following formula: 

T(Q,m,q) Q ~°° ^T/-f(q,m) * T {qjimi} -f(Q, m 7 , q 7 ) . (36) 

7 

Here, V is the Feynman diagram under consideration, {Q} ({m,q}) is the collection of the large (small) 
parameters, and the sum goes over all subgraphs 7 of T with masses m 7 and external momenta q 7 , subject to 
certain conditions to be described below. T{ q rn y is an operator performing a Taylor expansion in {q, m} before 
any integration is carried out. The notation r/7 * T!{„ m i7 indicates that the subgraph 7 of T is replaced by 
its Taylor expansion which should be performed in all masses and external momenta of 7 that do not belong 
to the set {Q}- In particular, also those external momenta of 7 that appear to be integration momenta in T 
have to be considered as small. Only after the Taylor expansions have been carried out, the loop integrations 
are performed. In the following we will refer to the set {7} as hard subgraphs or simply subgraphs, to {r/7} as 
co-subgraphs. 

The conditions for the subgraphs 7 are different for the hard-mass and the large-momentum procedure^. For 
the large-momentum procedure, 7 must 



• contain all vertices where a large momentum enters or leaves the graph 

Actually they are very similar and it is certainly possible to merge them into one condition using a more abstract language. 
For our purpose, however, it is more convenient to distinguish the two procedures. 



• ue une-par licit; ineuuauie aiiw luemiiyiiig mese vemces. 



From these requirements it is clear that the hard subgraphs become massless integrals where the scales are 
given by the large momenta. In the simplest case of one large momentum one ends up with propagator-type 
integrals. The co-subgraph, on the other hand, may still contain small external momenta and masses. However, 
the resulting integrals are typically much simpler than the original one. 

In the case of hard-mass procedure, 7 has to 

• contain all propagators carrying a large mass 

• be one-particle irreducible in its connected parts after contracting the heavy lines. 

Here, the hard subgraphs reduce to tadpole integrals with the large masses setting the scales. The co-subgraphs 
are again simpler to evaluate than the initial diagram. 

The large-momentum and hard-mass procedure provide expansions for off-shell external momenta which are 
either large or small as compared to internal masses. Recently a procedure allowing the asymptotic expansion 
of Feynman integrals near threshold was suggested. Graphical prescriptions similar to those for the large- 
momentum and hard-mass procedure have been worked out and applied to the cross section of quark production 



at e + e colliders near threshold 43 1. The expansion parameter in this case is given by the velocity of the 
produced quarks. 



A method to expand on-shell Feynman diagrams was developed in [81, 82]. Two typical limits for a large 
external momentum on a mass shell were considered: one where the mass shell is itself large, and the other one 
where the mass shell is zero. The latter case, called the Sudakov limit, is a purely Minkowskian phenomenon. 
This distinguishes it from the cases described so far which can be formulated completely in Euclidean space, 
simplifying rigorous proofs. The "philosophy" of these expansions is very similar to the hard- mass and large- 
momentum procedure. However, apart from the criteria on the subgraphs, the way of performing the expansion 
of propagators is also different. For example, the expansion of lines carrying large masses and not belonging 
to a one-particle-irreducible component of the subgraph (so-called cut lines) looks as follows: 

1 



K Kk 2 + 2Qk 



(37) 

K=l 



k is an integration momentum and Q a large external momentum. The graphical representation of the procedure 
becomes much less transparent because of that. 



In [|32[] the two-loop master integral with two different masses and on-shell external momentum was considered: 
m <C M, Q 2 = M 2 . The first 19 terms of the expansion in m/M were evaluated. Similar computations have 
also been performed for the fermion propagator j83|. Two-loop vertex diagrams with external momenta p\ 
and P2 obeying the Sudakov limit pf = and (p± + P2) 2 —00 were examined in [84|, and an expansion in 
m 2 /(pi +P2) 2 was obtained. 

Let us finally emphasize that due to analyticity the obtained expansions often provide valuable information 
also in other regions of the parameter space. This knowledge was used to reconstruct the photon polarization 
function by combining the results of asymptotic expansions in different limits (see also Section |4.1[). 



2.2.2 Examples 

Let us first consider the one-loop contribution to the photon propagator shown on the l.h.s. of the diagrammatic 
equation in Fig. [|. Both fermion lines are supposed to carry the same mass, m, and q is the external momentum. 
The application of the large-momentum procedure leads to the subdiagrams shown on the r.h.s. of this equation. 
The first one represents a simple Taylor expansion w.r.t. m, thus leading to massless one-loop integrals which 
can be solved with the help of Eq. (|3|). Starting from a certain order in m 2 , this subdiagram develops infra-red 
poles which are absent in the original diagram. They are due to the massless denominators arising in the Taylor 
expansion. However, they cancel against twice the ultra-violet poles of the second subgraph. The factor of two 
arises because a symmetric subgraph should be considered as well. To make this cancellation more transparent 



Figure 5: Large-momentum procedure for the one-loop photon polarization function. 



we present the first four terms of an expansion in m 2 for the transverse part of the one-loop polarization 
function (cf. Section fO]): 
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(38) 



with = la(—q /fjr), l qm = ln(—q /m ) and /^ m = ln(// /m ). The first two lines correspond to the first 
subgraph of Fig. |I| and the terms in square brackets are due to the second subgraph. The fourth and fifth lines 
show the sum of all subgraphs which corresponds to the consistent expansion of the full one-loop diagrams. The 
remaining pole in 1/e is the ultra-violet divergency of the full diagram and is usually removed by renormalization 
(see below). On the other hand, all spurious poles cancel in the sum. 

Let us now analyze the diagram on the l.h.s. of Fig. || in the limit q 2 <C m 2 . In this case the hard- mass 
procedure applies. It leads to a trivial Taylor expansion, and one ends up with bubble integrals. The first few 
terms read: 
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(39) 



Note that the 1/e pole is the same as in Eq. (|38|). In the case of the photon propagator the pole is usually 
removed by requiring that the polarization function vanishes for q 2 = 0. Finally the one-loop polarization 
function in the two limiting cases reads: 
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(40) 



Consider now the the case of the double-bubble diagrams pictured in Fig. |(| They provide a gauge invariant 
subclass of all three-loop graphs contributing to the photon polarization function. Note, however, that the 
particle type is irrelevant for the procedures described in the previous section; only the mass and momentum 
distribution is important. The outer and the inner mass are denoted by mi and ni2, respectively, and q is the 
external momentum. 

Of special interest is the high energy expansion of the current correlators, i.e., q 2 3> m^m 2 ., for example the 
case where mi = and mi = m. The imaginary part leads to 0(a 2 ) corrections to R(s) in the energy range 
sufficiently large as compared to the mass of the produced quarks. One can think of the production of light 



(a) (b) 
Figure 6: Fermionic double-bubble diagrams with generic masses mi and m2- 



quarks (u, d, s) accompanied by a pair of charm quarks (m = M c ) at energies yfs ;> 5 GeV, taking into account 
charm mass effects. 

Application of the large-momentum procedure to the diagram of Fig. |6](a) results in the subdiagrams displayed 
in Fig. |7[ We have combined topologically identical diagrams there. Similar terms arise from the graph of 
Fig. ^(b). Since the hard subgraphs have to be Taylor expanded in m and any "external" momentum except 
q, they belong to the class of massless two-point integrals, and can therefore be computed with the help of the 
integration-by-parts technique. The co-subgraphs, on the other hand, are massive tadpole integrals, so that 
also here integration-by-parts can be applied. Note that the first of the four subdiagrams in Fig. ^ represents 
the naive Taylor expansion of the integrand with respect to m. It is clear that starting from a certain order this 
term contains infra-red poles which are artificial as the original diagram is infra-red finite. The cancellation of 
these poles in the sum of all terms on the r.h.s. of the equation in Fig. ^ provides a non-trivial check of the 
calculation. The result of this expansion up to 0(m 8 /q 8 ) looks as follows: 
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(41) 



with I 

computation of all subdiagrams can be done using program packages designed to automate the large-momentum 



qui 



ln(-<f/W) and C3 = 1.202056903.... In Section [3.2.2| it is shown how the 



procedure. The logarithmic terms up to the fourth order can be found in 1 85] , the others are new. We have 
adopted the MS scheme, i.e., after taking into account the counter-terms for a s induced by the (massless) one 
and two-loop diagrams, the local poles are subtracted. Note that because the mass m is absent in the lower 
order terms, it does not need to be renormalized. 

As an application of the hard-mass procedure let us consider the hierarchy <C q 2 -C m 2 ,. The imaginary 
part again leads to contributions for R(s). This time one may think of charm quark production (mi = M c ) in 
the presence of a virtual bottom quark (1712 = Mj,). It turns out that already the first term provides a very 
good approximation almost up to the threshold = 2M5 [8C-88J. For simplicity we set mi = and 777-2 = 
in the following. 

The corresponding diagrammatic representation is shown in Fig. [8|. There are three subdiagrams, one of which 
again corresponds to the naive Taylor expansion of the integrand in the external momentum q. After Taylor 
expansion, the subdiagrams are reduced to tadpole integrals with mass scale m. The scale of the co-subgraphs 



is given by q thus leading to massless propagator-type integrals. The result for the first three terms reads [89|: 




Figure 7: Large-momentum procedure for the double-bubble diagram. The inner 
loop (thick lines) carries mass m, the outer one is massless, and so are the gluon 
lines. The square of the momentum q flowing through the diagram is supposed to 
be much larger than m 2 . It is understood that the hard subgraphs (right of "*") are 
to be expanded in m and all external momenta except for q, and reinserted into the 
fat vertex dots of the co-subgraphs (left of "*" ) . Contributions involving massless 
tadpoles are not displayed since they are zero in dimensional regularization. 
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Figure 8: Hard-mass procedure for the double-bubble diagram. Now the hierarchy 
q 2 <C m 2 is considered where m is the mass of the inner line. The hard subdiagrams 
(right of "*" ) are to be expanded in all external momenta including q and reinserted 
into the fat vertex dots of the co-subgraphs (left of ) . 
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If mi was different from zero, successive application of the large-momentum procedure to the co-subgraphs 
would lead to a subsequent expansion in m\jq 2 . An example for such a repeated use of asymptotic expansions 
will be described in Section |4.2j . Alternatively, one may evaluate the integration analytically as it has been 
done for the q 2 /m 2 term in |PQ| , thus leading to the full mi dependence of the power-suppressed result. 



2.3 Helicity-amplitude technique 

The standard way to obtain a cross section or a decay rate in perturbative quantum field theories is to compute 
the squared amplitude and integrate over the phase space for the final state particles. A typical amplitude is 
a sum over terms of the form 

c n e ^ (*.•> xi) n «(p; > A i ) r i u '(p'j - x j ) - ( 4s ) 

where e are polarization vectors of vector particles with momentum fcj and polarization x%i u an d u' are spinors 
of fermions (or anti-fermions) with momentum pj, p'- and helicity Xj, X'j, respectively. Tj are matrix- valued 
objects in Dirac space, and c is a scalar function of momenta, masses, coupling constants etc. Each term of the 
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square the amplitude before evaluating the expressions in ( f43| ) any further. After summing over polarizations 
of final and initial states, one employs the relations 

2_, u ijPi tyuip, A) ~ i> ± m (44) 

A=±l/2 

for fermions/anti-fermions and takes the trace in Dirac space which can be evaluated with the help of the anti- 
commutator {7^,71/} = iQpuv Specific polarization configurations can be investigated by introducing suitable 
projectors in Dirac space. 

However, the higher the order of perturbation theory, the more diagrams contribute, and because the amplitude 
must be squared, the number of terms to be evaluated in the above way increases even quadratically with the 
number of diagrams. As long as one stays with integrated quantities (total rates), one way out is to apply the 



optical theorem, i.e., to take the imaginary part of the forward scattering amplitude (see, e.g., Sections 4.1 



and 4.2). But as soon as one is interested in differential distributions, this method is no longer applicable. A 



solution here is the so-called helicity-amplitude technique. The basic idea is to rewrite expressions of the form 

U (p,X)u( P ',X') (45) 

for fixed helicities A, A' in terms of Dirac matrices. This allows the trace to be taken before squaring the 
amplitude. For given four-momenta one arrives at a single complex number for the full amplitude, one for each 
helicity configuration. If desired, the summation over polarizations is done only after squaring the amplitude. 



This strategy goes back to [91] and was further developed in [92]. Improved algorithms have been worked out 



in [93,94]. For detailed discussions of the various methods let us refer to these original works. 



2.4 Algebraic Programs 

Automatic computation of Feynman diagrams would not be possible without the development of powerful 
algebraic programming systems. While these systems themselves are based on conventional languages like 
C or LISP, in turn they constitute the breeding ground for all the application software which is our main 
concern in this review. Therefore, before discussing concrete packages that often are exclusively designed for 
Feynman diagram calculations, it may be helpful to briefly introduce four of the most frequently used algebraic 
programming languages, with emphasis on their different strategies and philosophy. 



2.4.1 Mathematica and MAPLE 

Mathematica certainly is the program with the largest bandwidth of possible applications. Apart from the 
algebraic operations, Mathematica performs analytical and numerical integrations for quite huge classes of 
integrals, it can do finite and infinite sums, differentiation, matrix operations, it provides graphic routines etc. 

MAPLE is very similar to Mathematica and often it is only a question of taste which one to use. There certainly 
are fields where one of them is superior to the other, and in extreme cases one may use them in combination 
to find the result of a calculation. 

The advantage of writing applications in Mathematica or MAPLE is that a lot of algorithms are built in and, of 
course, the whole environment — the interactive platform constituting an important piece in this respect - 
is available for the user to operate on the input and output. For example, one may perform numerical studies 
by inserting numbers for the parameters in the symbolic result. Meanwhile, a lot of users provide additional 
packages for very different purposes. 

One of the disadvantages of this generality is, of course, the rather big requirements of disk and memory space. 
Furthermore, the improper use of some very general commands (Integrate [] , Simplify [] , etc.) may lead to 
inadequately long operational times or even hang-ups, a problem mainly occurring with unexperienced users. 



2.4.2 FORM 



The power of FORM lies in its capability of manipulating huge algebraic expressions, the only limit in size being 
essentially the amount of free disk space. While for Mathematica or MAPLE operating on expressions of a few 



mega-uyies uaumes lamer inconveiiieiii, rurtn may easily ueai wim several nunuieu mega-uyies anu even giga- 
bytes. However, it is very restricted in its available operations. But after some experience one will be surprised 
how powerful the so-called id-statement is which replaces one quantity by an in general more complicated 
expression. 

The operational philosophy of FORM is completely different to the ones of Mathematica and MAPLE, and so 
far FORM is clearly optimized for algebraic operations in high energy physics. The concept is mainly to bring 
any expression to a unique form by fully expanding it into monomials ("terms"). Consider, for example, the 
following screen-shot of a Mathematica-session: 

In[l] := f 1 = (a+b)~2 

2 

Out[l]= (a + b) 

In [2] := f2 = a~2 + 2*a*b + b~2 

2 2 
Out [2]= a + 2 a b + b 

In [3] := diff = f 1 - f2 

2 2 2 

Out [3]= -a - 2 a b - b + (a + b) 

In [4] := Expand [diff] 

Out [4]= 

This clearly demonstrates that the internal representations of f 1 and f 2 inside Mathematica are not the same. 
FORM, in contrast, will instantly expand the bracket in f 1, i.e., the input file 

Symbol a,b; 
NWrite statistics; 

L fl = (a + b)~2; 

L f2 = a~2 + 2*a*b + b~2; 

L diff = fl - f2; 

print ; 
. end 

leads to the following output: 
fl = 

2*a*b + a~2 + b~2; 
f2 = 

2*a*b + a~2 + b~2; 
diff = 0; 

All further operations on f 1 or f2 will be term-wise, i.e., an id-statement 



id a = b + c; 
. sort 
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. sort means to collect all terms that only differ in their numerical coefficient and to bring the new expression 
for f 1 to a unique form again. This term-wise operation is the reason for FORM's capability of dealing with such 
huge expressions: Only a single term is treated in each operational step, all the others can be put on disk or 
wherever FORM considers it convenient. 

In addition, identification of certain terms ("pattern matching") is more transparent than in Mathematica. 
For example, the user must be aware of the fact that in order to nullify the mixed term in f 1, the function 
Expand [] has to be applied before. To be concrete, in Mathematica: 

In [2] := f 1 /. a*b -> 
2 

Out [2]= (a + b) 

In[3]:= Expand[ fl ] /. a*b -> 

2 2 
Out [3]= a + b 

It should be clear that the discussion above is not meant to judge if FORM's or Mathematica's concept is the 
better one. Each one of them has its pros and cons, although the number of operations in FORM is much more 
restricted than in Mathematica. 



2.4.3 REDUCE 

REDUCE [|ll[] is an algebra program that also was created for particle physics. It is not so general like 
Mathematica, but for specific problems especially in high energy physics it may be more efficient. In con- 
trast to FORM it has an interactive mode and many more built-in functions, but again its capability of dealing 
with huge expressions is more limited. 



3 Implementation 



In Section |3,l| a brief survey of the program packages is given which are very frequently used in high energy 
physics. Section |3.2| describes those packages in more detail which were used to obtain the physical results 
presented in Section ||[ 



3.1 Survey of the existing program packages 

This section gives a short overview of the existing program packages written to automate the treatment of 
Feynman diagrams. Meanwhile quite a number of such packages exist, some of them published in journals and 
available via anonymous ftp, others still under development and therefore only accessible for a restricted group. 
This review is not supposed to serve as a catalogue to select the one package most suitable for ones purposes. 
It shall simply provide an indication of what has been done, what is doable, and what are the most urgent 
topics to improve. Therefore we will not respect the question of availability. Furthermore, the description of 
packages and applications unavoidably will be biased, with emphasis on the topics the authors are and were 
personally involved in. However, since it allows us to go into details as far as ever necessary, we hope the reader 
will benefit from this strategy. Nevertheless, we will try to be as objective as possible and to fairly cover the 
class of most successful packages. 

The first part will be devoted to programs concerned with the generation of diagrams. For higher order correc- 
tions this becomes more and more important, because the number of diagrams increases rapidly with the order 
of perturbation theory. The second part describes programs that apply to the evaluation of the corresponding 
amplitudes on a diagram- by-diagram basis. Some of them are optimized for the use in combination with one of 
the generators mentioned before. Part three deals with two packages that automatically apply the rules of the 
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mainly combine some of the previously mentioned ones in order to treat full processes from the very generation 
up to the summation of the results of all diagrams. 

The aim of this section is not to compare or judge the listed programs. This is inadmissible anyway since each 
of them has its specific main focuses. We will just describe their needs, abilities and applications. Each package 
is based on one (sometimes also several) algebraic programs, the most important of which were introduced in 



Section 2.4. A summary in table form will also include some additional packages that could not be described 



in more detail. 



3.1.1 Generation of Feynman diagrams 



Two programs of quite different design will be discribed: FeynArts and QGRAF. 



FeynArts 

FeynArts [95] is written in Mathematica and may be used interactively. Given the number of loops and external 



particles, FeynArts first creates all possible topologies, allowing for additional criteria to select some subset 
of diagrams. In the second step, fields must be attributed to the lines which is most conveniently done by 
choosing a specific field theoretical model. The most popular models are predefined, but the user may equally 
well provide his own ones. FeynArts has the nice feature of drawing the generated diagrams which works to 
two loops with the default setup but may be extended to higher loop order by the user. This is very helpful 
for debugging, for example, by figuring out if the desired subset of diagrams was selected correctly. Finally, a 
mathematical expression is generated for each diagram which may then be further evaluated. This is preferably 
done using the program packages FeynCalc, FormCalc or TwoCalc (see below) because the user's intervention 
remains minimal in that case. It is very convenient that FeynArts works within the Mathematica environment 
since a lot of powerful commands are available here, allowing the manipulation of intermediate and final results. 



QGRAF 

The program QGRAF [^] is written in FORTRAN 77 and is a rather efficient generator for Feynman diagrams. It 
takes, for example, only a few seconds to generate 10000 diagrams. The user has to provide two files: the first 
one, called model file, contains the vertices and propagators in a purely topological notation. In the second 
one, the process file, the initial and final states as well as the number of loops must be defined. Furthermore, 
as for FeynArts, one may provide several options allowing the selection of certain subclasses of diagrams. 

Similar to the input, also the output of QGRAF is very abstract. It encodes the diagrams in a symbolic notation, 
thereby reproducing all necessary combinatoric factors as well as minus signs induced by fermion loops. While 
QGRAF suggests a distribution for the loop- and external momenta among the propagators, it is the users 
task to insert the Feynman rules, i.e., the proper mathematical expressions for the vertices and propagators. 
Nevertheless, QGRAF allows one to choose among different output formats which increases the flexibility for 
further operations. 



3.1.2 Computation of diagrams 



This section lists programs to be used for the computation of Feynman diagrams, all working on a diagram- 



by-diagram basis. Packages dealing with whole processes will be described in Section 3.1.4 



We will start with three packages based on FORM [13| and dealing with single-scale integrals: massive integrals 
with zero external momentum, massless integrals with one external momentum and two-point functions on 
their mass shell. Concerning their input requirements and their principle structure they are quite comparable. 

Furthermore we will describe certain Mathematica packages computing one- and two-loop diagrams and a 
MAPLE program which provides a graphical interface for such calculations. 



MINCER 



i_ne pacKage n±iM^Eiiv cuinpuies une-, iwu- aim unee-iuup imegiais wneie an lines are inassiess aim uniy une 
external momentum is different from zero. The first version of MINCER [p7| was written in SCHOONSCHIP 01, 



but here we will only describe the much more elaborate FORM- version [38|. MINCER was the first implementation 
of the integration-by-parts algorithm and certainly is one of the most important programs for multi-loop 
calculations. It is supposed to be highly efficient, in particular because the author of FORM was involved in the 
translation to this system. 

The user has to provide the diagrams that may all be listed in a single file, separated with the help of the fold- 
option of FORM. The input notation is based on the following scheme: The momenta carried by the propagators 
are numerated by integer numbers fixed by the topology of the diagram. The topology itself, uniquely classifying 
the diagram, must be specified through a keyword. The final result is given as an expansion in e where one-loop 
results are expanded up to 0(e 2 ), two-loop ones up to 0(e) and three-loop ones up to the finite part. 



MATAD 

One-, two- and three-loop vacuum integrals can be evaluated with the help of MATAD 8£], also written in FORM. 



Each propagator may either be massless or carry the common mass, M, and all external momenta have to be 
zero. 

The input notation as well as the whole concept is very similar to the one of MINCER. In particular, the core of 
the routines is again formed by the integration-by-parts algorithm. MATAD provides an interface for using MINCER 
and, furthermore, it supports Taylor expansion in small masses or momenta. For example, given a massive 
diagram with one small external momentum, after Taylor expansion of the integrand in the small momentum 
one is left with a vacuum integral again. At different stages of the computation there is the opportunity to 
interact from outside in order to control the expansion, apply projectors, or perform other operations. 

The applicability of both MINCER and MATAD seems to be quite restricted. However, as we have seen in Sec- 
tion if a diagram has a certain hierarchy of mass scales which happens to be the case for quite a lot of 
physical applications, it can be reduced to products of single-scale integrals. Then MINCER and MATAD can be 
used in combination to arrive at an analytical result. 



SHELL2 

The program SHELL2 |9^| is also written in FORM and deals with one- and two-loop propagator-type on-shell 
integrals. The implemented topologies allow computations mainly in QED and QCD. The prototype examples 
are the two-loop contribution to g — 2 of the electron and the relation between the MS and on-shell mass up 
to order a 2 in QCD. 

The user has to provide the diagrams in terms of a polynomial representing the numerator and abbreviated 
denominators. This must be supplemented by a label indicating both the number of loops and the topology. 
Except for the on-shell momentum no other external momentum may be present because then the numerator 
can be completely decomposed in terms of the denominators. A small FORM program ensures that an input 
very similar to the one of MATAD and MINCER can be applied so that one may run all three packages in parallel. 



FeynCalc, FormCalc and TwoCalc 



The Mathematica package FeynCalc [100] is based on a quite different philosophy than the ones described 
above. Its main applications are one-loop radiative corrections in the Standard Model and its extensions. 
The diagrams may either be provided by hand or one uses the output of the generator FeynArts (see above). 
FeynCalc performs the Dirac algebra and applies the tensor reduction algorithm of Section |2.1.3| in order 
to express the result in terms of scalar integrals. Special functions and the power of Mathematica allow to 
conveniently handle the intermediate and final expressions. 

Since FeynCalc is fully based on Mathematica, its performance is rather slow if the underlying expressions get 



large. Combining the advantages of Mathematica and FORM, the package FormCalc [ 101 1 is a sped-up version 
of FeynCalc well suited for huge problems at the one-loop level. 

The Mathematica package TwoCalc [61] is to some degree the extension of FeynCalc to two loops. It applies, 
however, only to two-point functions, since only for them may the tensor reduction algorithm be generalized 
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like the evaluation of Dirac traces, are passed to FeynCalc, so that mainly the reduction of the tensor integrals 
to a basic set of one- and two-loop integrals is performed. 

The numerical evaluation of the scalar one-loop integrals may conveniently be performed with the help of 
LoopTools to be described below. The two-loop integrals resulting from the TwoCalc-routines may be evaluated 
using the C-programs s21se and master | ]102| |. 



LoopTools pOl is an integration of the FORTRAN program FF [ |103[ into Mathematica. It allows a convenient 
numerical evaluation of the results as obtained by FeynCalc or FormCalc. In addition, it extends the ability of 
FF from doing only scalar integrals to the coefficients of the tensor decomposition By, Cy^, -Dyfc/ as described 
in Section 2.1.5 , so that it is not even necessary to fully reduce the tensor integrals with the help of FeynCalc 
in order to arrive at numerical results. Indeed, FormCalc reduces the tensor integrals only up to the point 
where LoopTools is applicable. 



ProcessDiagram 

The package ProcessDiagram [104] is written in Mathematica and deals with one- and two-loop vacuum 
diagrams without any restrictions concerning masses. One may apply Taylor expansions with respect to external 
momenta before integration. In addition, some auxiliary functions to simplify the results are available. It is, 
for example, possible to expand the result with respect to the masses. 



XLOOPS 



One 
XLOOPS 



and two-loop diagrams can be processed in a convenient way with the help of the program package 
An Xwindows interface based on Tcl/Tk allows even unexperienced users to compute loop 



lOg . 



diagrams. After choosing the topology, the particle type of each line must be specified. XLOOPS then performs 
the Dirac algebra and expresses the result in terms of certain one- and two- loop integrals. Whereas the one- loop 
integrals are evaluated analytically, for the two-loop ones the result is reduced to an at most two-fold integral 
representation which is then evaluated numerically. The main underlying strategy in this reduction is based on 
the parallel space technique [106 . For the analytical part of the calculations the algebra program MAPLE V [107 
is used; the numerical integrations are done with the help of VEGAS [108, 109]. 



3.1.3 Generation of Subdiagrams 

This section describes two packages, LMP and EXP, that automatically apply the hard-mass and the large- 
momentum procedure to a given diagram. The basic concept and the realization of these programs is quite 
similar. However, whereas LMP is written in PERL and was specialized for the large-momentum procedure, 
EXP is written in FORTRAN 90 and applies to the case of large external momenta, large internal masses, and a 
combination of both. 



LMP and EXP 



LMP [11C] is written in PERL and was developed to fully exploit the power of the large-momentum procedure. It 
therefore extends the range of analytically calculable three-loop diagrams from single-scale ones to two-point 
functions carrying small masses. It found its main application in the evaluation of quark mass corrections to 
certain QCD processes. 

Its basic concept was then carried over to the FORTRAN 90 program EXP [9C] which not only evaluates the 



large-momentum as well as the hard-mass procedure but also any successive application of both. Although the 
first version of EXP is only capable of dealing with two-point functions, the next version will also incorporate the 
case of an arbitrary number of external momenta and masses, so that any three-loop integral will be calculable 
in terms of a (possibly nested) series in ratios of the involved mass scales. So far, EXP was applied to the 



decay rate of the Z boson into quarks 25], but one can certainly think of a huge number of further possible 
applications. 

The usage of both LMP and EXP is very similar. Their input and output is adapted to MINCER and MATAD, so 
that the experienced user of these two packages will only have to provide some additional input information. 



In this section five program packages are described which automate the evaluation of a given process from the 
generation up to the computation of the corresponding amplitudes. Both their methods and purposes are quite 
different. 

GEFICOM computes Feynman diagrams up to three loops in analytical form by expanding them in terms of single- 
scale integrals; processes containing different mass scales are evaluated in terms of expansions. A different setup 
was used at NIKHEF (Amsterdam) to compute, for example, four-loop tadpole diagrams up to their simple 
poles in e. Further on, CompHEP is a multi-leg system calculating cross sections at tree level, involving up to five 
particles in the final state. The automation of GRACE relies heavily on numerical methods and has a similar field 
of application like CompHEP. Finally, an approach to provide a rather general environment for the computation 
of one- and two-loop diagrams is described with the programs TLAMM and DIANA. 



GEFICOM 

The program package GEFICOM [111| combines the generator QGRAF, the integration packages MATAD and MINCER 
and the programs EXP and LMP concerned with asymptotic expansion to compute Feynman diagrams up to three 
loops. The translation of the QGRAF output to MINCER/MATAD notation is done by a collection of Mathematica 
routines. The links between the single packages is done by script languages like AWK and PERL. 

With the short descriptions of its components in the previous section it is rather clear what the purpose of 
GEFICOM is: Given the initial and final states, GEFICOM generates and computes all contributing diagrams up 
to three loops in analytic form, provided that the involved mass scales are subject to a certain hierarchy. The 
result is obtained in terms of an in general multiple expansion in the ratios of the mass scales. The input 
therefore is essentially the same as for QGRAF, except that one additionally assigns masses to the different 
particles and defines a hierarchy among them. 



"NIKHEF setup" 

When the number of diagrams contributing to a single quantity increases and gets of the order 10 4 one should 
think carefully about the organization of the calculation. For example, storing the result of each individual 
diagram to a separate file may push the operational system over its limits. In the "NIKHEF setup" a database- 



like tool named Minos [112] is used to circumvent such book-keeping problems. It contains make-like and lots 
of additional features. For example, it helps to find bottlenecks of the setup by reporting on the subproblem 
on which most of the CPU time was spent. 

For the generation of the diagrams the program QGRAF is used. The output of QGRAF is translated to the notation 



of the FORM routines MINCER and BUBBLES [113] which are concerned with the integration of the massless two- 



point functions and the massive tadpole integrals, respectively. Another FORM program, named Color [114], 



determines the colour factor. The resulting expressions are inserted into a database. After integration, the 
results of the single diagrams are written to another database. Finally, the diagrams are multiplied by their 
colour factor and summed up. 



CompHEP 

CompHEP J!?]] is a program package which allows the evaluation of scattering processes and decay rates at tree 
level. A menu-driven interface makes its use quite handy. There are several built-in models among which one 
finds the Standard Model both in unitary and in 't Hooft-Feynman gauge, or the Minimal Supersymmetric 
Standard Model. Modification of these models and definition of new ones can be done in a very convenient 
way. Using LanHEP (see Section p.1.5 ) which works out the Feynman rules in CompHEP format, it even suffices 
to provide a Lagrangian density. 

After a model is selected and the process is specified, the Feynman diagrams are generated and graphically 
displayed on the screen. The user may now select the diagrams to be treated further. Then the squared 
Feynman amplitudes are generated and displayed, and the corresponding analytical expressions are computed. 
They may be stored either in REDUCE or Mathematica format which simplifies further symbolical manipulations. 
For complicated processes FORTRAN or C code is generated, allowing for numerical studies. On the other hand, 
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showing angular distributions and cross sections can be produced. 



The numerical part of CompHEP is based on the Monte-Carlo integration routine VEGAS. It is possible to introduce 
cuts on various kinematical variables. Further on, distributions, cross sections and particle widths can be 
evaluated. For the incoming particles one may define structure functions and then repeat the same integrations. 
Finally, CompHEP generates events and displays the corresponding histograms. 



GRACE 



A different approach for the automatic computation of Feynman diagrams is realized by the program package 
GRACE [115]. It was developed to compute cross sections and radiative corrections to them. Also here several 
models are available, including the Minimal Supersymmetric Standard Model [116]. 

GRACE has its own generator for Feynman diagrams, producing FORTRAN source code for each of them and 
passing it to the package CHANEL []117| which performs the calculation with purely numerical methods. The 



user is free to directly evaluate the amplitudes rather than the squared ones. This significantly reduces the size 
of the expressions (see Section ^^ ). 

In the final step the integration over the phase space for the particular final state is performed with the help 
of the multi-dimensional integration package BASES and the event generator SPRING which allows to generate 
unweighted event flow [115]. 



TLAMM and DIANA 



There are two projects, partly under development, the purpose of which is to automate the evaluation of 
processes involving one- and two- loop diagrams. The package TLAMM [US] is written in C and was mainly 
developed to compute the two-loop corrections to the anomalous magnetic moment of the muon. It uses QGRAF 
to generate the diagrams, translates the output to FORM code and executes the integration routines. 

The second package, called DIANA [)120f | , is supposed to work on a more universal basis. Again, in a first step 
the output of QGRAF is read, the topologies are determined and internal representations for the diagrams are 
created. Subsequently, an interpreter executes a "special text manipulating language" (TM) which allows one 
to select the algebra language to be adopted for the calculation, to pass the expressions to FORTRAN in order to 
perform a numerical calculation or to generate, for example, PostScript files of the diagrams. 



3.1.5 Miscellaneous 



There are many more programs developed by several different groups. Many of them are neither published 
nor documented in the literature, especially if they were designed for a very special task only. This section is 
supposed to touch on some of them in a more or less encyclopedic formP]. For completeness also the programs 
already discussed are included in the list below. Note that when we did not include information on the 
availability of a program this means that we could not find an apropriate ftp or http site. Interested readers 
should contact the authors of the corresponding programs in this case. Some programs are also available from 
the CPC program library (http://www.cpc.cs.qub.ac.uk) upon request, even if it is not indicated below. 



BASES ITT 



— Availability: CPC program library 

— Purpose: Monte Carlo integration 

— Source: FORTRAN 



BUBBLES JTT3J: 



— Purpose: analytical computation of purely massive four- loop tadpole integrals up to C(l/e) 

— Algorithms: integration- by-parts 

— Source: FORM 



2 The links to the corresponding www and ftp sites can also be found at 
http : //www-ttp . physik.uni-karlsruhe . de/Links/ algprog.html 



uniiiM ill . 



Availability: CPC program library 

Purpose: library for the calculation of helicity amplitudes 
Source: FORTRAN 



Color [114]: 



— Availability: http: //norma. nikhef .nl/~t68/F0RMapplications/Color 

— Purpose: computation of colour factors 

— Source: FORM 



CompHEP J17): 

— Availability: http://theory.npi.msu.su/~comphepor http://www.ifh.de/~pukhov 

— Purpose: symbolic and numerical computation of tree level processes with up to six external legs 

— Algorithms: symbolic evaluation of squared diagrams, recursive representation of kinematics, Monte 
Carlo integration 

— Source: FORTRAN, C 

— Uses: VEGAS 

— Preferably combined with: LanHEP 
DIANA 1T2 



— Availability: upon request from the author 

— Purpose: general environment for generating and evaluating Feynman diagrams 

— Source: C 

— Uses: QGRAF 



EXP gO I 



— Purpose: reduce two-point functions to single-scale integrals 

— Algorithms: large-momentum procedure, hard-mass procedure (see Section [2.2| ) 

- Preferably combined with: MAT AD, MINCER 

- Source: FORTRAN 90 



FeynArts Q: 

— Availability: http : //www-itp . physik . uni-karlsruhe . de/f eynarts 

— Purpose: diagram generator with main focus on one- and two-loop cases 

— Source: Mathematica 

— Preferably combined with: FeynCalc, FormCalc, TwoCalc 



FeynCalc flOO[ : 



— Availability: http://www.mertig.com (commercial); 

the latest free version is available at http://www.mertig.com/oldfc 

— Purpose: reduction of arbitrary one-loop to a set of basis integrals 

— Algorithms: tensor decomposition and tensor reduction by means of |14j (see Section 2.1.3) 

— Source: Mathematica 

— Preferably combined with: FeynArts, LoopTools 



FF ||64|103|: 

— Availability: http://www.xs4all.nl/~gjvo/FF.html 

— Purpose: numerical computation of scalar and vector one-loop integrals up to six-point functions 

— Source: FORTRAN 



FormCalc [101 : 



Availability: http : //www-itp . physik . uni-karlsruhe . de/f ormcalc 

Purpose: slimmed, high speed version of FeynCalc 
Algorithms: tensor decomposition 
Source: Mathematica, FORM 

Preferably combined with: FeynArts, LoopTools 



- Purpose: automatic generation and computation of three-loop Feynman diagrams in terms of ex- 
pansions 

- Uses: QGRAF, MATAD, MINCER, EXP or LMP 

- Source: Mathematica, FORM, AWK 

GRACE ill 



— Availability: ftp://ftp.kek.jp/kek/minami/grace 

— Purpose: numerical computation of 2 — > 2 scattering processes to one-loop and multi-particle scat- 
tering processes at tree level 

— Algorithms: helicity-amplitude method (see Section |2.3| ), Monte Carlo integration 

— Uses: CHANEL, BASES, SPRING 

— Source: C, FORTRAN 



HELAS [121] : 



Purpose: helicity-amplitude subroutines for Feynman diagram evaluations 
Source: FORTRAN 



HEPLoops J12J: 

— Availability: upon request from the author 

— Purpose: analytical computation of massless propagator-type diagrams up to three loops 

— Algorithms: integration-by-parts (see Section 2.1.2| ) 

— Source: FORM 



LanHEP [p3fl : 

— Availability: http://theory.npi.msu.su/~semenov/lanhep.html 

— Purpose: generate Feynman rules from Lagrangian 

— Preferably combined with: CompHEP 

— Source: C 



LMP |110| : 

— Purpose: factorize large external momentum in two-point functions 

— Algorithms: large-momentum procedure (see Section |2.2| ) 

— Source: PERL 

— Preferably combined with: MATAD, MINCER 



LOOPS |L2J: 

— Availability: CPC Program Library 

— Purpose: computation of one- and two-loop propagator type integrals 

— Algorithms: integration-by-parts (see Section |2.1,2j ) 

— Source: REDUCE 



LoopTools []101|] : 

— Availability: http : / / www-itp . physik . uni-karlsruhe . de/looptools 

— Purpose: implementation and extension of FF in Mathematica 

— Uses: FF 

— Source: Mathematica, FORTRAN 

— Preferably combined with: FeynCalc, FormCalc 



MadGraph |125([ : 

— Availability: http : //pheno . physics . wise . edu/Sof tware/MadGraph/ 

— Purpose: automatic generation of Feynman diagrams; calculation of helicity amplitudes 

- Uses: HELAS 

- Source: FORTRAN 



master [102|: supplement to s21se 



— Purpose: analytical computation of massive three-loop tadpole integrals 



— Algorithms: integration- by-parts (see Section 2.1.2 ) 

- Source: FORM 

MINCER (FORM version) 



— Availability: ftp : //nikhef h.nikhef . nl/pub/theory/f orm/libraries/f orm2/mincer 

— Purpose: analytical computation of massless propagator-type diagrams up to three loops 

— Algorithms: integration-by-parts (see Section 2.1.2) 

— Source: FORM 

MINCER (SCHOONSCHIP version) 

— Availability: CPC program library 

— Purpose: analytical computation of massless propagator-type diagrams up to three loops 

— Algorithms: integration-by-parts (see Section 2.1.2) 

— Source: SCHOONSCHIP 



MINOS fll2fl : 

— Purpose: controlling facility for the calculation of processes with a huge number of diagrams 

— Source: C 



oneloop [126|: 

— Availability: http://wwwthep.physik.uni-mainz.de/~xloops 

— Purpose: algebraic and numerical calculation of one-loop diagrams 

— Source: MAPLE 



PHACT fl27H : 

— Purpose: numerical computation of tree processes up to four particles in the final state 

— Algorithms: helicity amplitudes by means of |)4| 

— Source: FORTRAN 



ProcessDiagram [104]: 

— Purpose: computation of one- and two-loop bubble diagrams allowing for several different masses 

— Source: Mathematica 



PVEGAS |p.09|l : parallel version of VEGAS 

— Availability: ftp://ftpthep.physik.uni-mainz.de/pub/pvegas/ 

— Source: C 



RECURSOR p28[ : 

— Purpose: analytical computation of massive three-loop tadpole integrals 

— Algorithms: integration- by-parts (see Section 2.1.2) 

— Source: REDUCE 

s21se (10 

— Availability: f tp : //f tp . physik . uni-wuerzburg . de/pub/hep/index . html 

— Purpose: numerically evaluate integral representations for scalar two-loop self-energy integrals 

— Source: C 

— Preferably combined with: TwoCalc 
SHELL2 Hf: 

— Availability: CPC Program Library 

— Purpose: computation of propagator-type on-shell integrals up to two loops 

— Algorithms: integration-by-parts (see Section 2.1.2) 

— Source: FORM 



SIXPHACT ||129||: extension of PHACT to six particles in the final state 



— Purpose: analytical computation of massless propagator-type diagrams up to three loops 

— Algorithms: integration- by-parts (see Section 2.1. 2| ) 

— Source: REDUCE 

SPRING HUgJ : 

— Purpose: event generation 

— Source: FORTRAN 



TARCER plf : 

— Availability: http://www.mertig.com/tarcer/ 

— Purpose: tensor reduction of two-loop propagator-type integrals and reduction to a set of basic 
integrals 

— Algorithms: tensor reduction by shifting the space-time dimension (see Section |2.1,4 ) 

— Source: Mathematica 



TLAMM |TTgl: 

— Purpose: automatic evaluation of two-loop vertex diagrams 

— Source: C 



TRACER |L3| : 



Purpose: Dirac trace calculations in Mathematica, optionally with 't Hooft-Veltman or naive anti- 
commuting 75 

Source: Mathematica 



• TwoCalc Kl: 



— Availability: upon request from the author 

— Purpose: reduction of two-loop propagator diagrams to a set of basis integrals 

— Algorithms: two-loop tensor reduction by means of [61] 

— Source: Mathematica 

— Uses: FeynCalc 

— Preferably combined with: FeynArts, s21se, master 



• QGRAF |6[: 

— Availability: ftp://gtae2.ist.utl.pt/pub/qgraf/ 

— Purpose: efficiently generate multi-loop Feynman diagrams in symbolic notation 

— Source: FORTRAN 

• VEGAS fllTT 



— Availability: see, e.g. [ 133 

— Purpose: adaptive multi-dimensional integration 

— Source: FORTRAN 



WPHACT [|3J : extension of PHACT to W and Higgs physics 

— Availability: http://www.to.infn.it/~ballestr/wphact/ 



XLOOPS p05| : 

— Availability: http://wino.physik.uni-mainz.de/~xloops/ 

— Purpose: compute one- and two-loop diagrams using X-interface 

— Algorithms: parallel space technique [106] 

— Uses: PVEGAS, oneloop 

— Source: MAPLE, Tcl/Tk 



We will start with the description of the FORM programs MAT AD and MINCER. Details on the automatic application 
of asymptotic expansions using LMP and EXP are given in Section 3.2,2 . In Section 3.2.3 the package GEFICOM 
will be explained. The way the generator QGRAF works should also become clear then, which is why it will not 
be described separately. The subsequent section is devoted to the Mathematica packages FeynArts, FeynCalc, 
FormCalc, TwoCalc and LoopTools. In the concluding section we discuss the package CompHEP as an application 
of automated multi-leg computations. 



3.2.1 MATAD and MINCER 

Let us examine the three-loop ladder diagram (see Fig. ^) in order to see in detail how the programs MATAD 
and MINCER work. We will focus on the case where the external momentum, q\, is much smaller than the mass 
of the fermions, m. Application of the rules for the hard-mass procedure (Section [T^) shows that it suffices to 
perform a Taylor expansion of the integrand in order to arrive at an expansion in q 2 /m 2 . The resulting integrals 
to be evaluated are of the tadpole type, and the calculation can completely be passed to MATAD. In the opposite 
limit, q\ 3> m 2 , in addition to the naive Taylor expansion of the integrand in m 2 /q 2 , the large- momentum 
procedure generates non-trivial subgraphs. In general their evaluation requires to use MATAD and MINCER in 
combination. 

P2 




Figure 9: Ladder-type diagram contributing to the photon propagator. The solid, 
wavy and curly lines represent quarks, photons and gluons, respectively. 



The input for the Feynman graph is written to a file using the fold construct of FORM 



* — #[ ladder: 
((-1) 

*Dg(nul,nu2,p7) 
*Dg(nu3,nu4,p8) 

*S(mul ,ql ,p3m,nu3 ,ql ,p2m,nul ,ql ,plm,niu2,p6ni,nu2 ,p5m,nu4,p4m) 
); 

#define TOPOLOGY "LA" 
* — #] ladder: 

Some words concerning the notation are in order: The function S ( . . . ) represents the fermion string and 
contains momenta (pi =pi, ...) and Lorentz indices (mul = fit, ■■■) as its arguments. The momenta 
are labeled in accordance with the MINCER conventions for the specific topology which is LA here, meaning 
"ladder" . Note that it is more convenient to express the input in terms of eight momenta and to use momentum 
conservation at a later stage than to impose it already in the input. The ending m indicates that massive fermions 
should be considered. A small FORM program transforms each Lorentz index to the corresponding 7 matrix 
and the momenta to fermion propagators. The expansion in the external momentum, ql, is indicated by the 
consecutive appearance of ql together with an integration momentum. The depth of the expansion, which is 
encoded as a pre-processor variable in FORM, must be defined at program call. Dg(. . .) represents the gluon 
propagator. After reading the diagram these functions are replaced according to the Feynman rules. Then 
MATAD evaluates the fermion trace and expands the denominators with respect to ql. At this point the user 
may influence the procedure by providing an additional file, the so-called "special treat-file". For the matter 
in hand, this file contains a projector on the transverse part of the diagram: 



13|: 



iuu.j. L±pj.y , ^ _ <-i_ \mu.i ,mu.zj -rq± ^mm,/ fqi \.muzj / q± . q±,/ ^ueno , -z.) ; 



where deno(3,-2) means 1/(3 — 2e). In this way one obtains a scalar expression where ql is absent in the 
denominators of the propagators. However, it may still appear in scalar products with integration momenta. 
Application of d'Alembertian operators performs an averaging in the numerator, i.e. terms like (q\ • p) 2 are 
replaced by p 2 q 2 /D. This projects out the coefficients of (ql) 71 which have no dependence on the external 
momenta any more. After that it is possible to decompose the numerators in terms of the denominators. At 
this stage the expression may already have split into millions of different integrals. Recurrence relations derived 
with the integration-by-parts technique reduce them to a small set of tabulated master integrals which have been 
evaluated once and for all. In the case of the photon propagator, for example, all three-loop integrals reduce 
to only three massive three-loop master integrals. The same is true for the massless case with non-vanishing 
external momentum, which is where MINCER applies. 

In intermediate steps, when the recurrence relations are applied, it may happen that the expression blows up 
and disk space up to several Gigabytes is needed. But the final result is rather compact and fills only of the 
order of one output screen (depending, of course, on the depth of the expansion). 



3.2.2 EXP and LMP 



The programs EXP and LMP are computer-technical realizations of the prescriptions for the asymptotic expansion 



of Feynman diagrams described in Section 2.2. They both work at the graphical level, i.e. they are not dealing 



with the analytic expressions corresponding to the Feynman diagrams but rather with the diagrams themselves. 

Given an arbitrary diagram with a single external momentum (not counting momenta that can be gotten rid 
of by a simple Taylor expansion), EXP and LMP generate all possible subdiagrams by successively removing any 
combination of lines (propagators). For those subdiagrams matching the conditions of the large- momentum, 
respectively, the hard-mass procedure, the programs determine the momentum distribution. In particular, 
they determine which lines should be Taylor expanded in which momenta, and how the momenta of the hard 
subgraph are related to the ones of the co-subgraph. Thus, every diagram is split into several subdiagrams, all 



of them being products of single-scale integrals as discussed in Section ]2J2 



As it was already mentioned in Section p. 1.3 , LMP is written in PERL and is only concerned with the large- 



momentum procedure. EXP, being in some sense its successor, is a FORTRAN 90 program^ and deals also 
with the hard-mass procedure. In addition, if several scales are involved, it repeatedly applies both methods in 



combination, thus reducing an arbitrary three-loop self-energy diagram to single-scale integrals (see Section 1.2). 



To be specific, let us consider the example for the double-bubble diagram of Section 2.2.2 where the outer line 
is massless and the inner one carries mass m. The input file for LMP looks as follows: 

*--#[ GLOBAL : 
#define POWER "12" 
*~#] GLOBAL : 



*--#[ TREAT : 

multiply, (-d_(mu,nu) + q(mu)*q(nu)/q.q)*deno(3,-2) ; 
*~#] TREAT : 



*--#[ db : 

*S(mu,pl , ro ,p2 ,nu,p3 , si ,p4) 
*SS (tau , p6m , al , -p7m) 
*Dg(tau,ro ,p8) 
*Dg(si,al,p5) 

#define TOPOLOGY "01" 
*--#] db : 



where the fold construct of FORM is used not only to encode the diagram itself but also the projectors (TREAT- 
fold) and possible options (GLOBAL- fold). The projector is again on the transverse part of the diagram (see 

3 A forthcoming version will be written in C. 



Decuuii |o.z.i| ;. iiie tJAuerxiai liiumeiiiuiii is uenuieu uy q mis nine. i_iie uniy upuuii nere is ruwun, ueiiuuiig 
the required depth for the expansion in the mass m. 

The diagram is contained in the fold named db (for "double bubble"). The encoding is based on MINCER/MATAD 
notation (see Section 3.2.1| ). In particular, the momentum distribution is in accordance with the corrsponding 



topology defined in these programs. 

LMP reads the content of the db-fold, translates it to an internal graph-based notation and applies the large- 
momentum procedure. There are two main output files; the first one is called db.dia and contains the relevant 
subdiagrams: 



*--#[ db_l : 

*S (mu , +p 1 , r o , +p2 , nu , +p3 , s i , +p4) 
*SS (tau , +p6mexp , al , -p7mexp) 
*Dg(tau,ro,+p5) 
*Dg(si,al,+p5) 

> 

#define DIANUM "1" 
#define TOPOLOGY "arb" 
#define IMT1 "inpol" 
#define MASS1 "0" 
*--#] db_l : 



[. . .] 



*--#[ db_5 : 

*S(mu,+pll,ro,+pl2,nu,+pl3,si,+pl4) 
*SS (tau,plm, al ,+aexp, -pl5mexp) 
*Dg(tau,ro,+pl5) 
*Dg(si,al,+pl5) 

#define DIANUM "5" 
#define TOPOLOGY "arb" 
#define IMT1 "topLl" 
#define MASS1 "M" 
#define INT2 "inptl" 
#define MASS2 "0" 
*--#] db_5 : 



[. . .] 



Only two out of six contributing subdiagrams are displayed here, the remaining ones being represented by 
[ . . . ] . The information on the topology of the diagrams is now split according to the factorization of loops. It 
is contained in the pre-processor variables INT1, INT2, etc. The variable TOPOLOGY is set to the dummy value 
arb. 

The first subdiagram, db_l, corresponds to the naive Taylor expansion of the integrand with respect to m. This 
can be seen from the fact that the momentum distribution is the same as in the original diagram, and that the 
massive lines are denoted by pimexp instead of pirn. In this way MATAD is advised to perform a Taylor expansion 
in m. The second subdiagram displayed above, db_5, corresponds to the second one on the right-hand side 
of Fig. [?] (note that there is a topologically identical one, indicated by the factor 2 in Fig. [?]; LMP, however, 
generates both of them as separate diagrams). The co-subgraph is a one- loop tadpole-diagram carrying the 
momentum p\. It must not be expanded in the mass, therefore its propagator is simply denoted by plm above. 
The hard subgraph is a two-loop diagram of the topology shown in Fig. |]. Its momenta pi (i = 1, ... ,5) 
are denoted by pli, where the additional "1" is introduced to distinguish them from the momenta (here it is 
only one momentum, actually) of the co-subgraph. One propagator carries the momentum pi — pi§ which is 
indicated by the notation +aexp,-pl5mexp. Note that this propagator also has to be expanded w.r.t. m. Here, 
pi is denoted by aexp to let MATAD know that it has to expand with respect to p\. After expansion, aexp is 



luenuneu wiun pi. ±ms is une ui ine main cuiuems ui uie secunu umpui ine ui uir, caneu LieaL.au. ±i again 
splits into folds, one for each subdiagram. The fold corresponding to the subdiagram db_5 above, for example, 
contains the line 



id aexp = + pi; 

which does the identification of momenta mentioned above. Most of the remaining statements of the treat . db- 
file are concerned with the proper "power-counting", i.e., to take care that the Taylor expansions are performed 
such that not too few and not too many terms are generated. (Expanding up to unnecessarily high power may 
drastically slow down the computation.) 



To compute the second example of Section |2.2,2| , the one concerned with the hard-mass procedure, one has 
to use EXP. The input file, however, is almost identical to the one shown above. The only difference is the 
GLOBAL- fold, since in the case of EXP also the hierarchy of scales has to be fixed (for LMP, only q 2 S> m\ , m 2 ., . . . 
is allowed): 

*--#[ GLOBAL : 
#define POWERM "6" 
#define SCALES "M,q" 
*~#] GLOBAL : 

This means to assume m 2 3> q 2 and to expand up to (q 2 /m 2 ) 3 . For multi-scale problems, this fold may look 
like 

*--#[ GLOBAL : 
#define POWERMa "4" 
#define POWERMc "2" 
#define POWERQ "8" 
#define SCALES "Ma,q,Mc,Mb" 
*~#] GLOBAL : 

indicating the hierarchy m 2 S> q 2 S> m 2 ^> m 2 and an expansion up to the corresponding degrees. The output 
files of EXP are very similar to the ones above, although the power-counting statements are obviously much 
more involved than for LMP. 

In addition, LMP and EXP produce the whole set of administrative files, namely FORM files which control the 
running of MINCER/MATAD, and make files concerned with program calls and the reduction of loss caused by 
system crashes. 

The only task to be done by the user is to feed in the input file shown above, call LMP or EXP, and finally to 
type "make" in order to start the very computation of the diagrams. The output is exactly a result like the 
one quoted in Eq. ( |4l| ) and (^), respectively. 



3.2.3 GEFICOM 

This section investigates in more detail how to use GEFICOM to evaluate a whole set of Feynman diagrams. 
Fig. [n] gives a schematic overview of GEFICOM and the possible interactions from outside. It unifies several 
packages described above and thus allows to automatically compute a specified class of diagrams in analytic 
form as long as a certain hierarchy among the involved mass scales is justified. The strategy then is to 
apply asymptotic expansions in order to split all integrals into massless propagator-type diagrams and massive 
tadpoles. 

As an example let us consider corrections to the gluon propagator from which follows, for example, the gluon 
wave function renormalization constant. GEFICOM needs a name for each problem, and we will choose "gg" 
for it. A more sophisticated example involving up to four different dimensional parameters is considered in 



Section 4.2 



In a first step the user sets up a file containing the Lagrangian, or rather its topological content, i.e., the plain 
vertices and propagators without the corresponding analytical expressions. In our example it is called gg . lag 
and looks as follows: 
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Figure 10: The structure of GEFICOM. 



* propagators 

[g,g;3+t] 

[c,C;l-] 
[q,Q;2-] 

[sigma, sigma; 1+t] 

* vertices 

[Q.q.g] 
[C.c.g] 

[g.g.g] 
[g,g, sigma] 

This file serves as input for QGRAF^J and the notation is adopted from this program. In the first part the 
propagators are defined for gluons (g), ghosts (c) and quarks (q). The auxiliary particle a (sigma) is introduced 
to split the four-gluon vertex into three-linear vertices which are much easier to deal with. In addition to the spin 
multiplicity one also has to specify if the particle is a boson ("+") or fermion ("—")• The optional parameter 
t guarantees that no diagrams with tadpoles of the corresponding particle are generated. The notation for the 
vertices in the second part is rather obvious. 

A second file, called gg.def , contains the underlying process, accompanied by additional options: 



*** MINCER 
* gauge 1 



4 The notation of the QGRAF version 1.3 has been adopted. 



* power 

list = symbolic ; 
lagfile = ' q.lag' ; 
in = g [q] ; 
out = g [q] ; 
nloop = ; 
options = ; 

true = bridge[ g,c,q, 0,0 ]; 

The first line determines that the output of the generation procedure should be transformed into MINCER 
notation. The next two statements fix the gauge and the regularization scheme, and the variable power defines 
the depth of the expansion in small masses and momenta. Since in our case all particles are taken to be 
massless, only one scale is left, namely the external momentum. Thus, an expansion is not needed and the 
variable power may be assigned the value zero. The remaining part again uses QGRAF notation and specifies 
the diagrams to be generated. For example, "in" and "out" denote in- and outgoing states, and the last line 
excludes diagrams that are one-particle reducible w.r.t. gluons, ghosts or quarks, options is left empty, and 
nloop, the number of loops, will be defined via a command line option later. 



Like in Section 3.2.1, the user again may provide "special treat-files" — in the case of the gluon propagator 
two of them. The first one, called treat . gg. 1, again contains the projector on the transverse structure of the 
correlator: 

1 / _ q^q^ 2 



D-U^Ti' (46) 

Besides the consideration in momentum space GEFICOM also takes care about the colour factor which is why in 
general one should supply a projector also in colour space. For the gluon propagator with the colour indices 
a(l) and a(2) it is given by 

ra(l)a(2) 

, (47) 

n g 

where n g = N% — 1 is the number of gluons. The corresponding FORM file, called treatcol .gg. 1, looks as 
follows: 

multiply, prop(a(l) ,a(2))/ng; 

After providing these simple input files, GEFICOM may be called by executing the shell-script command: 
> doq2f -p <prb> -1 <loops> 

where <prb> gives the name of the problem (gg here) and <loops> is the number of loops. GEFICOM calls QGRAF 
and transforms the output to Mathematica format. Unfortunately QGRAF does not provide the topology of the 
generated diagrams. While for a human being it is very simple to figure out this information by just looking 
at the diagram, for a computer this is a non-trivial task, especially if the number of loops exceeds two. One of 
the core parts of GEFICOM indeed is a Mathematica program precisely concerned with this problem. 

Once the topology is available, the notation for the diagrams is translated to a format which is suitable for 
the FORM packages MINCER and MATAD. In addition, a set of administrative files controlling the calculation and 
ensuring minimal loss in the case of computer breakdowns is generated. For the case of the gluon propagator 
considered in this section, at three-loop level 494 diagrams are generated, and it takes roughly 12 hours of 
CPU time on a DEC-Alpha machine with 600 MHz for the complete evaluation. A sample diagram for the 



three-loop gluon propagator is shown in Fig. 11 



In analogy to the momentum-space diagram also a FORM version of the "colour diagram" is generated which is 
evaluated with the help of a small FORM program. 



Figure 11: Sample diagram contributing to the gluon propagator at three loops. 



The computation of the diagrams is not initiated automatically. Instead, GEFICOM generates a shell-script, 
do<prb> (dogg in our case) that allows, for example, only a subset of diagrams, or even only a single diagram, 
to be computed via command line options. The advantage of this strategy is that different computers may be 
used in parallel for different subsets, or that a part of the diagrams may even be evaluated before the generation 
is finished. After the results of all diagrams are available they are summed up, taking into account the correct 
colour (and other) factors and the (unrenormalized) result is at hand. 



3.2.4 FeynArts, FeynCalc, FormCalc, TwoCalc and LoopTools 

FeynArts: FeynArts is a Mathematica-package generating Feynman diagrams and the corresponding ampli- 
tudes. Its great flexibility has allowed it to become an enormously useful tool not only for high energy particle 
physicists. It has found its applications also in effective field theories, and by knowing its power one could 
imagine that it may be useful even for at first sight completely unrelated subjects. 

There exists a manual describing its features in a rather clear and pedagogical way. The four following functions 
constitute the heart of FeynArts: CreateTopologies, InsertFields, Paint and CreateFeynAmp. 

As an example let us consider the one-loop corrections to the triple Higgs vertex in the Standard Model. The 
use of CreateTopologies in this case is the following: 

tops = CreateTopologies [1 ,3, Tadpoles -> False, Self Energies -> False]; 

The first two arguments of CreateTopologies determine the number of loops and external legs, respectively. 
The remaining arguments are options, where the first one prevents creation of tadpole insertions and the second 
one states that no diagrams with self-energy insertions on external legs are generated^. To view the topologies, 
one types 



Paint [tops ,1,2] ; 



where the first argument specifies the list of topologies and the second and third one refer to the number of 



incoming and outgoing fields, respectively. The outcome is the graphic shown in Fig. 12 
Next the fields are attributed to the lines. This is done by typing 



ins = InsertFields [tops, {S[l]> -> {S [1] , S [1] } , Model -> {SM}] ; 



where the first argument again refers to the list of topologies created above, the third one fixes the model (SM 
= Standard Model), and the second one defines the in- and outgoing particles (S[l] for the physical Higgs 
boson). This produces a list of 64 diagrams that can be viewed using Paint again, this time, however, with 
only one argument, since in- and outgoing fields are already specified in "ins": 



Paint [ins] ; 



An extract of the list of diagrams is shown in Fig. 

5 In what follows we will use the notation of the FeynArts version 1.2 

6 There is a single option comprising both of the latter, named WFCorrections (for "wave function corrections"). 
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Figure 12: Topologies created by FeynArts. 
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Figure 13: Sample of diagrams generated by FeynArts for the triple Higgs vertex 
to one-loop order. 



The last step to be done by FeynArts is to generate analytic expressions which are understood by FeynCalc, 
FormCalc and TwoCalc: 

amps = CreateFeynAmp [ins] ; 

This produces a list containing the amplitudes of the individual diagrams. The 17 th element corresponds to 
the diagram with a top-quark triangle (Fig. |14|) and reads as follows: 

In [12] := amps [[17]] 

3 1 3 3 

0ut[12]= FeynAmp[SlSlSl, Tl, 117, N17] [(— EL MT Integral [ql] 

128 

> tr[(MT + gs[kl - pi + ql]) . (MT + gs [ql] ) . (MT + gs [kl + ql])]) / 

34222 22 23 

> (MW Pi (-MT + ql ) (-MT + (kl + ql) ) (-MT + (kl - pi + ql) ) SW )] 

pi are the ingoing, kz the outgoing, and qi the loop momenta. EL is the electric charge, gs [p] means ^ and tr 
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Figure 14: 17 th dia gram produced by FeynArts for the triple Higgs vertex. 



the Dirac trace. 



FeynCalc, FormCalc and TwoCalc: The list of expressions named amps above may now be directly fed into 
FeynCalc]] or FormCalc. Its entries will automatically be transformed to the internal notation (see below). 
The key function of these programs is OneLoop, respectively derivatives of it, in particular OneLoopSum. As a 
simple example consider the integral 

/ (L)V?-? 2 )K-(^) 2 ) = is? ■ ^ 

In the notation of FeynCalc this reads 

In [3] := OneLoop [q, l/(2*Pi) ~4 FeynAmpDenominator [ 
PropagatorDenominator [q,ml] , 
PropagatorDenominator [q+k,m2] ] ] 

I 2 2 

— BO[k.k, ml , m2 ] 
16 

Out [3]= 

2 

Pi 

Note that the pre-factor in In [3] above is ^-independent. In principle, this could be assigned to a definition 
of BO in Out [3] that differs from the one in (El) by a D-dependent factor; however, another interpretation is 
to implicitely assume MS-regularization of all integrals. In fact, as will become clear in a moment, this is what 
is done in the numerical routine LoopTools. 

In this way FeynCalc rewrites any one-loop diagram to the standard integrals defined in Eq. (|9|). For example, 
the diagram corresponding to amps [[17]] defined above produces the output 

In [4] := OneLoopSum [amps, SelectGraphs -> {17}] 
Out [4]= K[7] 

In [5] := FixedPoint[ReleaseHold [#]&,%] 



7 We use version 2.2/3 here which is the latest version which is freely available (see Section 



3.1.5). 



-3 EL MT BO[kl.kl, MT , MT ] 3 EL MT BO [pi. pi, MT , MT ] 

Out [5]= 

3 2 3 3 2 3 

32 MW Pi SW 32 MW Pi SW 

3 4 2 2 

3 EL MT BO[kl.kl - 2 kl.pl + pl.pl, MT , MT ] 



3 2 3 
32 MW Pi SW 

3 4 2 2 2 

(3 EL MT CO[kl.kl, pl.pl, kl . kl - 2 kl . pi + pi . pi , MT , MT , MT ] 

2 3 2 3 

(4 MT - kl.kl + kl.pl - pl.pl)) / (32 MW Pi SW ) 

The action of FormCalc is quite similar to the one of FeynCalc, except that is does not perform the full 
reduction to the scalar integrals of Eq. (|Tl|), but rather only decomposes the expressions into covariants, 
keeping the corresponding coefficients as they are. In general this reduces the length of the expressions and 
still can be evaluated numerically with the help of LoopTools, to be described below. The same strategy can 
also be followed by FeynCalc, simply by setting the option ReduceToScalars in OneLoop to False. The main 
feature of FormCalc in comparison to FeynCalc is its speed which is mainly due to the fact that it passes large 
expressions to FORM. 

TwoCalc is the extension of FeynCalc to two-loop propagator type diagrams, and its usage is very similar. 



FF and LoopTools: FF is a FORTRAN program which numerically evaluates the scalar one-loop integrals defined 
in ([lT]). The extension to the coefficients appearing in the tensor decomposition of tensor integrals and its 
integration into Mathematica is done in LoopTools. E.g., the numerical evaluation of the i?o-mtegral defined 
above may be simply done in Mathematica by saying 

In[2]:= B0[k.k,ml~2,m2~2] //. {k.k -> ml~2,ml -> 100, m2 -> 10} 
Out [2]= -7.49109 

This example shows that the 1/e poles (together with the constants 7e and ln4-7r) are set to zero as it corresponds 
to an MS-renormalized quantity; this was already pointed out in the discussion below Eq. (|48|). Both the ultra- 
violet and the infra-red divergent parts are controlled with the help of parameters which specify the square 
of the renormalization scale /U. In a similar way, it is possible to obtain numerical values for the diagram 
amps [[17]] above by providing numbers for the scalar products and masses. 



3.2.5 CompHEP 



Let us now discuss in more detail how the package CompHEP works. As already mentioned in Section |3.1.4 



CompHEP essentially consists of two parts: the symbolic and the numerical one. A flowchart representing the 



symbolic part is shown in Fig. |15 



After choosing the underlying model in menu 1, one gets to menu 2 where the process must be defined and 
certain diagrams may be discarded from further considerations with the help of menu 4- One may also modify 
the parameters of the selected model (see menu 3). menu 5 is concerned with the graphical output and the 
algebraic squaring of the matrix elements for the generated diagrams. 

In menu 6, analytical expressions for the squared matrix elements are computed. They can be stored to disk 
either in C, FORTRAN, REDUCE or Mathematica format (see menu 8). Instead of applying the helicity-amplitude 
technique, CompHEP uses the conventional method of squaring the amplitudes before calculating traces. Whether 
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Figure 15: Menus of the symbolic part of CompHEP [1351 
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Figure 16: Main menu of the numerical part of CompHEP [135] 



this is an advantage or a disadvantage depends on the specific problem under consideration. For simple processes 
a built-in numerical interpreter (see menu 7) is available to perform the phase space integration. For more 
complicated problems one has to initiate the actual numerical part of CompHEP. It deals with processes that 
exceed the capability of the built-in calculator. 



The main options of the numerical part of CompHEP are displayed in Fig. |16[ Before starting the computation, 
the environment for the Monte Carlo integration (using VEGAS) must be set. This includes, for example, the 
modification of model parameters and the choice of the center-of-mass energy. Furthermore, the user decides 
whether structure functions for the incoming particles should be used or not. A central point is represented by 
the item Regularization, where the integration variables, previously defined in the item Kinematics, are mapped 
in such a way that the integrand becomes a smooth function. This is one of the main features of CompHEP: 
The users task is just to provide possible singularities in the propagators that are close to the edges of the 
phase space region; the mappings themselves are done completely automatically. Other items are concerned 
with cuts over energy, angles, transverse momenta, squared momentum transfers, and invariant masses or the 
rapidity for a set of outgoing particles. Note that all the numerical calculations may equally well be performed 
in batch mode which is useful as soon as the run-time becomes large. 



In order to reduce the CPU time, CompHEP saves the matrix elements computed during a VEGAS run into a 
file. This file is short because it contains only one number for each event. Later-on, a special program called 



gencveiiLB repeals me iviome ^ano evaiuauoiis, leaunig uie previously caicuiaieu maunx elements iruni uie 
file instead of performing a re-calculation. Thus the CPU time is drastically reduced for this second run, being 
typically of the order of only a few minutes as compared to a few hours for the initial run. Therefore, it is 
possible to study the influence of various additional cuts on the cross-section and to fill histograms rather 
quickly by performing several subsequent runs of genEvents. Each working session produces two output files: 
one containing the results of the calculation together with a list of model parameters and a copy of the screen 
report, and a second one containing the computed matrix elements. 



4 Applications 

This section reports on recent results that would have been impossible to obtain without the use of computer 
algebra. The selected examples are supposed to underline the fields of applications for some of the programs 
described above. 

The first example is concerned with one of the most classical subjects for multi-loop calculations, the photon 
polarization function. We will describe its evaluation, using the program MAT AD in the limit q 2 <C m 2 , and 
additionally MINCER and LMP in the opposite case, q 2 S> m 2 . The technique of repeatedly applying the hard- 
mass procedure was successfully applied to the decay rate of the Z boson into b quarks using EXP in the 
GEFICOM environment. The calculation of renormalization group functions was the driving item for BUBBLES 
to be developed. As part of the "NIKHEF-setup" it was used to compute four-loop tadpole diagrams up 
to their simple poles in e. Another example for the calculation of a four-loop quantity is the decay rate of 
the Higgs boson into gluons. Applying the approach of an effective Hamiltonian, GEFICOM could be used to 
compute the corresponding diagrams. As an application of the multi-leg program CompHEP the scenario of a 
strong interaction among electroweak gauge bosons will be discussed. It is traded as an alternative to the Higgs 
mechanism to restore unitarity at high energies. Calculations in the full electroweak theory usually involve a 
large number of different scales. For important contributions to the parameter Ar it was possible to reduce 
the contributing diagrams to two-point functions which then were accessible with the help of the programs 
FeynArts and TwoCalc. 



4.1 QCD corrections to the photon polarization function 



Notation and methods 

This section is concerned with the computation of QCD corrections to the photon propagator. We define the 
polarization function in the following way: 

+ n(g 2 ) = iJd 4 xe^(0\Tj^x)j u (0)\0), (49) 

where = ipj^ip is the diagonal vector current of two quarks with mass m. The main motivation is the simple 
connection of Ii{q 2 ) to the physical quantity R{s) which is defined as the normalized total cross section to the 
production of heavy quarks: 

R(s) ee " ( ff - ha f ons ) = 12 nImU( q 2 = s + ie). (50) 

The advantage of making a detour to H(q 2 ) and not considering R(s) from the very beginning is the implicit 
summation of real radiation and infra-red singularities. The fact that the number of loops to be evaluated for 
H(q 2 ) is larger by one is made up for by having to deal with two-point functions only, for which a huge amount 
of technology is available. 

It is convenient to separately consider the following contributions: 

ii(q 2 ) = n(°)(, 2 ) + ^c F n( 1 )(g 2 )+f^Vn( 2 )(g 2 ) + ... , (51) 

IT* 2 ) = C 2 F U ( ^ +C A C F U^ A + C F Tn l U < i 2) + C F TUP, (52) 



Figure 17: Diagrams contributing to the one-, two- and three-loop polarization 
function. Solid lines represent quarks, loopy ones are gluons. 



and similarly for R(s). The colour factors Cp = (N 2 — 1)/(2N C ) and Ca = N c correspond to the Casimir 
operators of the fundamental and adjoint representations of SU(3), respectively. The case of QCD corresponds 
to N c = 3, the trace normalization of the fundamental representation is T = 1/2, and the number of light 

I — K (2) 

(massless) quark flavours is denoted by n/. In Eq. (|52[), 11^ is the Abelian contribution (corresponding to 



quenched QED) and IT]^ is the non-Abelian part specific for QCD. There are two fermionic contributions 

(2) 

arising from diagrams with two closed fermion lines, so-called double-bubble diagrams: For the quark in 

(2) 

the inner loop is massless, the one in the outer loop massive, whereas for Up both fermions carry mass m. 
The case where the external current couples to massless quarks and these via gluons to massive ones will not 
be addressed here. Its imaginary part was considered in IB 7[|. The Feynman diagrams contributing to one-, 



two- and three-loop order are shown in Fig. 17 



The two- loop corrections of 0(a s ) were computed in analytic form in the context of QED quite some time 



ago [136]. In subsequent works the calculation was redone and more convenient representations were found [137 — 
|139| . It is, at least with the currently available technology, out of range to compute the three-loop diagrams 
in complete analytic form for arbitrary q and m. Only for a subclass of diagrams, namely the ones with two 
closed fermion lines where one of them is massless, was it possible to compute the corresponding contribution 
to R{s) analytically []87, 140| , 88, 141]. Nevertheless, there are essentially two approaches which provide very 



good approximations for U^ 2 \q 2 ), both of them relying heavily on the use of computer algebra. 

One method is to consider the polarization function n(g 2 ) in the limit q 2 3> m 2 . Application of the large- 
momentum procedure (see Section |2.2| ) leads to an expansion in m 2 /q 2 with the coefficients still depending on 
logarithms ln(— q 2 /m 2 ). The idea is to evaluate as many terms as possible (the limitations essentially coming 
from the CPU time and the size of the intermediate expressions) and thus to approximate the true function 
Tl(q 2 ) even for rather small values of q 2 . 

Whereas the above method results in analytical expressions, an alternative way is to construct a semi-numerical 
result for II (q 2 ) which will, however, be valid for all values of q and m. In addition to the high energy expansion 
discussed above, this method also requires the expansion of Tl(q 2 ) in the limit q 2 <C m 2 . Together with some 
information about the (two-particle) threshold q 2 = Am 2 , a suitable conformal mapping and the use of Pade 



approximation, one obtains an approximation for II(g 2 ) over the full kinematical region |2C, 142 1. 

Referring to the literature for physical applications (e.g., [ 143] , 144]), in the following we will concentrate on 
the technical aspect of the calculation in the limits q 2 <C m 2 and q 2 S> m 2 . While the latter case directly 
corresponds to the first method above, the former one is needed for the second method. 

As can be seen from Fig. |l7| only a small number of diagrams contribute at three-loop level. This makes the use 
of a diagram generating program like QGRAF unnecessary. The main challenge instead is to compute as many 
terms in the expansions as possible. Apart from fast computers this requires an efficient algebraic language 
and optimized programs to deal with a large number of terms. The choice for this problem was MATAD and 



ri-UM^Eitt, uum wimen m rurm i^see Decuuii p.4.i\ ). 



Large mass limit 

In the limit q 2 <C m 2 application of the hard-mass procedure (see Section |2.2| ) shows that it suffices to keep 
the naive Taylor expansion of the diagrams in their external momentum. Thus one stays with the calculation 
of massive tadpole integrals. Let us consider the moments C n of the polarization function, defined through the 
Taylor series 

n« 2 V) = ^c^UX, (M) 

n>0 \ / 

where m is the on-shell mass. Although this series does not develop an imaginary part, one may gain information 
on the rate from it by exploiting the analyticity of Il(q 2 ). We do not want to discuss this approach here and refer 



the interested reader to [145, 20|. At three-loop level the evaluation of the coefficients up to Cg was performed 



in [ 20 1 . The calculation required disk space of the order of several GB for FORM (so-called "formswap") to store 



the intermediate expressions. The total CPU time on a 256 MHz DEC-Alpha workstation with 128 MB main 



memory was roughly two weeks. With the input for the ladder-type diagram already shown in Section 3.2.1 
let us now have a look at the output of MATAD. The expansion up to terms of order (g 2 /m 2 ) 4 reads: 

ladder = 

+ ep~-3 * ( - 8/9*Q.Q ) 

+ ep~-2 * ( - l/693*Q.Q~5*M~-8 + 8/945*q . q~4*M~-6 - 2/35*q . q~3*M~-4 + 
8/15*q.q~2*M"-2 - 134/27*q.q ) 

+ ep"-l * ( - 20*M"2 - 4042699/108056025*q.q~5*M~-8 + 273004/1488375* 
q.q~4*M~-6 - 63029/66150*q.q~3*M~-4 + 10658/2025*q . q~2*M~-2 - 473/81* 

q.q - 4/3*q.q*z2 ) 

+ 96*z3*M~2 - 482/3*M"2 - 19508/225*q . q~5*z3*M"-8 - l/462*q . q~5*M~-8*z2 
+ 129728592122581/1248047088750*q.q~5*M~-8 + 42496/675*q . q~4*z3*M~-6 
+ 4/315*q.q~4*M~-6*z2 - 559517166977/7501410000*q . q~4*M~-6 - 44* 
q.q~3*z3*M~-4 - 3/35*q.q~3*M~-4*z2 + 6119442979/125023500*q.q~3*M~-4 

+ 328/27*q.q~2*z3*M~-2 + 4/5*q . q~2*M~-2*z2 - 880963/121500*q . q~2* 
M"-2 + 11909/486*q.q - 392/9*q.q*z3 - 67/9*q.q*z2; 

Note that this expression still contains an overall factor q 2 (the results above are in Euclidean space, indicated 
by the capital Q, i.e. Q 2 = —q 2 ). After dividing by q 2 there are terms of 0(m 2 /q 2 ) which cancel in the sum of all 
three-loop diagrams. The constant terms disappear after requiring the QED-like on-shell condition n(0) = 0, 
and finally the structure of Eq. ( |53|) is obtained. 

After all diagrams up to three loops are computed and added up (taking into account the correct colour factors) , 
the parameters a s and m are renormalized. To give an impression of the structure of the final result, we list 

(2) 

the first two and the eighth term of the Abelian part, 11^ : 



_C2) 8687 32 A n n A , 22781 

C\ { = C2 In 2 + 4 C2 H C3 

A ' x 864 5 s 1728 S 



223404289 192 t , 24 . 4857587 



°A,8 



190302182417255312898886115648452691 

63063833203636585050931200000 
786432 98304 31209476560803609727258477 

~ 230945 ^ 2 n + 46189 ^ + 12432176686773043200 ^ 3 ' 

where Cn = C( n )> with Riemann's ^-function. For the specific values above it is (2 = 7r 2 /6 and £3 = 1.20206 . 
The left-out coefficients and the ones for the other colour structures of (|5"2|) can be found in . 



CO * 



Figure 18: Sixth subdiagram as generated by LMP for the large-momentum proce- 
dure of the ladder diagram (cf. Fig. ||). 



Large momentum limit 



The calculation for q 2 3> m? is more involved in the sense that here the asymptotic expansion does not reduce 
to a naive Taylor expansion. Instead, the large-momentum procedure to be applied in this case generates a 



large number of subgraphs for each of the diagrams of Fig. 17. Whereas for the small-momentum expansion 



a computer only had to be used for the computation of the diagrams, in this case even the generation of the 
subdiagrams is a non-trivial task and program packages like EXP or LMP (see Section 3.2.2) are indispensable. 
For the three-loop case the 18 initial diagrams produce 240 subgraphs and their manual selection would be 
very tedious and error-prone if not impossible. The corresponding numbers in the two- (one-) loop case are 2 
(1) initial and 14 (3) subdiagrams which would still allow a manual treatment. These numbers also show that 
there is a big step in complexity when going from two to three loops. 



For the calculation the package LMP (see Section 3.2.2 ) was used. As an example let us again consider the 
ladder diagram of Fig. ^. The input file, looking similar to the one of the low-energy case above, produces 
27 subdiagrams. The one of Fig. 18 (the sixth subgraph generated by LMP), for example, corresponds to the 
following code: 



*--#[ d313_6 : 
((-1) 

*Dg (nul , nu2 , +aexp , -p2) 
*Dg (nu3 , nu4 , -bexp , +p2) 

* S (nu2 , +p lmexp , nu4 , -bexp , +qmexp , mu2 , -p2 lm , nu3 , -p2mexp , nu 1 , 

-pi lm , mul , -aexp , +qmexp) 
*D ; 

#define DIANUM "6" 
#define TOPOLOGY "arb" 
#define INT1 "inpll" 
#define MASS1 "0" 
#define INT2 "topLl" 
#define MASS2 "M" 
#define INT3 "topLl" 
#define MASS3 "M" 
*--#] d313_6 : 



After Taylor expansion, the integral factorizes into a massless and two massive one-loop integrals. The topology 
of each of these one-loop integrals is encoded in the pre-processor variables INT1, INT2, and INT3, where inp. . 
denotes a massless, top. . a massive topology (see also Section 3.2.2j ). The corresponding mass is given by 
MASS1, MASS2, and MASS3. 

(2) 

Again the first two and the highest available (in this case the seventh) term in the expansion for will be 
displayed: 
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(55) 



with l qm = ln(— q 2 /fh 2 ) and Z 9/J = ln(— q 2 //J. 2 ), ( n as in (|54|) and the additional values £ 4 = 7r 4 /90, ^5 = 
1.03693. ... -E>4 was analytically calculated in |5^] and evaluates numerically to B4 = —1.76280. ... m is the 
MS mass which can easily be transformed to the on-shell scheme |>f|. For the left-out terms in ( |55|) we refer 
to p6| , [l47| . 

In Fig. |l^ (a) the imaginary part of the Abelian contribution is plotted as a function of x = 2m/y/s including 
successively higher orders in x. Up to x ~ 0.7 very quick convergence is observed. For larger values of 
x, however, the inclusion of higher terms does not significantly improve the approximation. This becomes 
manifest in an even more drastic way in Fig. [l9| (b), showing the same curves as functions of v = y/1 — 4m 2 Js. 
This choice of variable enlarges the threshold region and thus demonstrates the breakdown of the high-energy 
expansion close to threshold. 
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4.2 Corrections of (D(otot s ) to the decay of the Z boson into bottom quarks 

In the examples considered so far essentially only the QCD part of the Standard Model is involved, which is, 
as already noted, the most important field for multi-loop calculations. The electroweak sector, on the contrary, 
resides on a broken gauge symmetry which is the reason why the diagrams in general carry a lot more scales. 
Higher order calculations are thus much more involved. 

As an in some sense intermediate case, one may consider mixed electroweak/QCD corrections. An important 
example in this class are the 0(aa s ) corrections to the decay rate of the Z boson into quarks. One may 
distinguish two cases, namely the decay into the light u, d, s or c quarks and the one into bb. The latter case is 
of special interest as the top quark enters the vertex diagrams already at one-loop order and leads to corrections 
proportional to M 2 . The complete 0(a) corrections were considered in | |148j,|149| |. Since the mixed 0(aa s ) 
correction to the Z decay into u, d, s and c turned out to be quite sizable |150(| , it was tempting to consider 
also the decay into bottom quarks at this order. The leading terms of 0{GfM 2 ) and 0(Gf ln(M f 2 /M^)) were 



computed in |151 - 154 



However, the full Mt dependence to 0(aa s ) is currently out of reach. On the other hand, the top quark is 
certainly the heaviest particle entering the Zbb vertex at this order, which makes asymptotic expansions a well 
suited tool to obtain a very good approximation to the full answer. 




(d) 





Figure 20: Diagrams containing a top quark that contribute to Z — > 66. Thin lines 
correspond to bottom quarks, thick lines to top quarks, dotted lines to Goldstone 
bosons and inner wavy lines represent W bosons. 



Let us describe the technical aspects of the calculation for Z — » bb in more detail. The 0(aa s ) corrections are 
computed by applying the optical theorem, i.e., by evaluating the imaginary part of the Z boson propagator. 
Our concerns are only the diagrams involving the top quark. The ones contributing to 0{a) are shown in 
Fig. [2(]. The diagrams to be evaluated to 0{aa s ) can be obtained by attaching a gluon in all possible ways, 
thereby increasing the number of loops by one. 

The different scales in these diagrams are given by Mt , A% , Mz and M$ . The latter one is the mass of 
the charged Goldstone boson. It is related to Myy via M| = ^w^w-> where £w is the corresponding gauge 
parameter which is arbitrary and should drop out in the final result. In a first step we will consider it such that 
M t 2 > M|. For the dia grams (a), (b), (d) and (e) the application of the hard- mass procedure with respect 
to Mt immediately leads to a complete factorization of the different scales, i.e., only single-scale integrals with 
at most three loops are left. For the diagrams (c), (f) and (g), however, some co-subgraphs still carry more 
than one scale which makes an evaluation quite painful, especially as their 0(e) part is also needed. It is very 
suggestive to apply the hard-mass procedure again to these co-subgraphs. In order to avoid unwanted imaginary 
parts arising from W and cuts the hierarchy to be chosen among the physical masses is Mf 3> M§. 

The last inequality is seemingly inadequate, but a closer look at the corresponding diagrams shows that in this 
context it is actually equivalent to 4M^ 3> M§, respectively, (Mjy + M t ) 2 3> M§. The mass of the Goldstone 
boson may be incorporated using the possible hierarchies £yyMyy S> M^ 3> M\ or ^> £ w Myy 3> M§. 
This procedure leads to a nested series in Myy/Mf and M z /M^r. Diagrammatically, an example for this 
procedure of successively applying asymptotic expansions looks as follows: 



+ 



+ 




where only those terms are displayed which are relevant in the discussion above and all others contributing to 
the hard-mass procedure are merged into the three dots. The mass hierarchy is assumed to be M t 2 £y(/M|L 3> 



The technical realization of the calculation was performed with the help of GEFICOM (see Section 3.2.3). The 
main new ingredient in comparison to the examples quoted so far is the use of the automated version of the 
hard-mass procedure and its successive application, implemented in the program EXP, in combination with the 
generator QGRAF and the FORM packages MINCER and MATAD. This automation not only avoids human errors, but 
also allows several checks of the results. For example, since the gauge parameter £w is arbitrary, it is possible 
to choose it such that no longer the top quark, but rather the Goldstone boson is the heaviest particle involved 
in the process. Then one is left with the hierarchy ^w^w ^ ^z- Note that for this choice 

the hard-mass procedure produces completely different subdiagrams than for the cases above, where Mt was 



eoiisiuereu 10 ub me laigesi acme, r ux i>xiurxu.um, wmie a manual ueaimem ui uum cases wuuiu aimusi uuuuie 
the effort, the automated version allows one to go from one hierarchy to the other by simply interchanging two 
input parameters, as we will see below. 

To be concrete, let us consider the input file needed for GEFICOM. It is certainly more sophisticated than the 
examples mentioned so far, but this is mainly due to the enriched particle spectrum. 



*** MINCER 



* scheme 2 

* gauge 1 

* exp y ma mc mb q 

* powerma 2 

* powermb 4 

* powermc 4 

* mass t Ma 

* mass Wp Mb 

* mass pp Mc 

* loops 1 true=iprop [t ,0,0] ; 

* loops 2 true=iprop [g,0,0] ; 

* loops 2 f alse=iprop [t ,0,0] ; 

* loops 3 true=iprop [g, 1 , 1] ; 

* loops 3 f alse=iprop [t ,0,0] ; 



list = symbolic ; 

lagfile = ' q.lag' ; 

in = Z [q] ; 

out = Z [q] ; 

nloop = ; 

options = onepi; 

true = iprop [Wm , pm , , 2] ; 

false = iprop [bq,0, 1] ; 



The meaning of most of the lines can be deduced by comparison with the example of Section |3.2,3| . The mass 
of the top quark (t =t), the W (Wp =W + ), and the Goldstone boson (pp =$ + ) are denoted by Ma, Mb and 
Mc, respectively. The change in mass hierarchy mentioned before is achieved by rewriting the line 



* exp y ma mc mb q 



to 



* exp y mc ma mb q 



This input file is combined together with the one containing the vertices and propagators, which shall not be 
displayed here, in order to generate in a first step the relevant diagrams at three-loop level. After inserting 
the Feynman rules, EXP applies the hard-mass procedure and reduces all diagrams to single-scale integrals, if 
necessary by applying the procedure twice. EXP also produces the relevant administrative files which then call 
MINCER and MATAD. The runtime for QGRAF is of the order of a few seconds and it takes a few minutes for EXP to 
generate the subdiagrams. The time spent in the integration routines strongly depends on the required depth 



of the expansion, of course. The computation in [25] took about three weeks. 

The contribution of the vertex corrections to the decay rate of the Z boson to bottom quarks, induced by the 
exchange of a W or a Goldstone boson, may be written in the following way: 

STV = SFf + {5Tf - 5Tf) = 5Tf + « W " K' W ) = STf + 6TZ d , (56) 



where 5Y^ and ^r^'^ denote renormalized and unrenormalized quantities. In the difference between the b and d 

b-d 



contributions the relevant counter-terms drop out. This means that ST]^ d is independent of the renormalization 



scneme aim is meieiuie wen suneu iui iiisiaiiauuii 111 uaw analysing piugiains line me unes uescnueu, e.g., 
155]. In addition, 5T^_ d is also gauge independent. 5T^ has been computed in [15C[. For convenience we 



in 



list the result for 5T^_ d in numerical form [25]: 
5Tf_ d = r ^ -{ - 0.50 + (0.71 - 0.48) + (0.08 - 0.29) + (-0.01 - 0.07) + (-0.007 - 0.006) 

Sn 7T I 



+ — 

IT 



Y — - 

2 

si vr 



1 a 



1.16 + (1.21 - 0.49) + (0.30 - 0.65) + (0.02 - 0.21 + 0.01) + (-0.01 - 0.04 + 0.004) 
1.16 + 0.13 



0.50 -0.07 + — 

7T 



(57) 



' and Mt the on-shell 
0.223. The numbers 



with T° = N c Mza/ (12sqCq) , sq = sin#vF> &w being the weak mixing angle, cf, = 1 — s 
top mass. For Eq. (|7|) we used the values M t = 175 GeV, M z = 91.19 GeV and sj 
after the first equality sign correspond to successively increasing orders in 1/M 2 , where the brackets collect 
the corresponding constant, ln(M 2 /M^) and, if present, In 2 (M 2 /M^) terms. The numbers after the second 
equality sign represent the leading M 2 term and the sum of the subleading ones. The O(a) and 0{aa s ) results 
are displayed separately. Comparison of this expansion of the one-loop terms to the exact result of []149| shows 
agreement up to 0.01% which gives quite some confidence in the aa s contribution. One can see that although 
the M t 2 , M t ° and M t ° ln(M 4 2 /M^) terms are of the same order of magnitude, the final result is surprisingly well 
represented by the M? term, since the subleading terms largely cancel among each other. 



4.3 Four-loop (3 function 

The (3 function is an object of common interest in any field theory, especially in non-Abelian ones. It describes 
the dependence of the corresponding coupling constant with respect to the energy scale. For QCD it is 
convenient to define 

^^IT = = -1^1 EM^I " (58) 

The (3 function is directly related to the renormalization factor Z g defined through = Z 2 a s , where and 
a s are the bare and renormalized coupling constant of QCD, respectively: 




P(a s ) = ^ 1 ^Z 2 / 1 l (59) 
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d 



it da s 



9 



where Zg'^ is the residue of Zg with respect to its Laurent expansion in e. In Eq. (^) it is already indicated 
that in mass independent renormalization schemes, like the MS scheme, the renormalization constants and thus 
also the (3 function only depend on a s . 

There are several ways to compute Z g , all of them related through Slavnov- Taylor identities. One choice would 
be to combine the renormalization constants of the four-gluon vertex and the gluon propagator. At four-loop 
level, however, this requires the computation of about half a million diagrams. On the other hand, the approach 
of [23] was based on the relation 

Zg = -J±=, (60) 

^3V Z 3 

where Z\ is the renormalization constant of the ghost-ghost-gluon vertex, and Z3 and Z3 are the ones of the 
ghost and gluon propagators, respectively. Thus the pole parts of the corresponding Green functions had to 
be computed up to four loops. 

The roughly 50,000 contributing diagrams were generated by means of QGRAF. Introducing an artificial mass M 
in each propagator as infra-red regulator allowed a Taylor expansion to be performed with respect to all external 
momenta, leading to four- loop massive tadpole integrals. Recurrence relations based on the integration- by-parts 



algorithm were implemented in the program BUBBLES [US] to reduce the integrals to a minimal set of master 



integrals. At one-, two- and three-loop level only one, at four-loop level two master integrals are required. 
The colour factors of the individual diagrams were determined with the help of a FORM program To 



cope wiiii me nuge nuiiiuei ui uiagrams a special uaiauase-iiKe iooi caneu i'iximud was ueveiupeu. ±i conuois 
the calculation and allows to conveniently access the results of single diagrams. The total CPU time for this 
calculation was of the order of a few months. 



The final result which we will only quote for the physically most interesting case of QCD reads 
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where nt is the number of active flavours. 

4.4 Hadronic Higgs decay 
4.4.1 Introduction 



As was already pointed out in the introduction, the only particle of the Standard Model not yet discovered 
is the Higgs boson. Thus it is necessary to learn as much as possible about potential production and decay 
mechanisms of such a particle in order to correctly interpret the signals in the detector. An important decay 
channel of the Higgs boson is the one into gluons. Although it is a loop-induced process it is numerically 
not negligible mainly due to the fact that the 0(a s ) corrections are very large [21,p2|. In this section it is 
shown how to compute the O(o^) corrections exploiting the powerful tools of automatic Feynman diagram 
computation. For completeness also the decay into quarks will be discussed. 

The Higgs boson will be assumed to be of intermediate-mass range, i.e., Mh & 2Myy, which suggests the use 
of the approach of an effective Lagrangian where the top quark is integrated out. For convenience of the reader 
let us, in a first step, collect the relevant formulae. Detailed derivations can be found in |i5^ , ^57| M[i58|. 



The starting point is the Yukawa Lagrange density describing the coupling of the H boson to quarks. In the 
limit — > oo it can be written as a sum over five operators [ 156j , 157] formed by the light degrees of freedom, 
with the corresponding coefficient functions containing the residual dependence on the top quark: 
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q£{u,d,s,c,b,t} 
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eff 



ttO 5 

i=l 



(62) 



It turns out that only two of the operators, in the following called 0[ and 0' 2 , contribute to physical processes. 
Expressed in terms of bare fields they read 



0[ = (G^ 



O' 



0/ 



q£{u,d,s,c,b} 



0/ 



(63) 



where G^ a is the gluonic field strength tensor. The renormalization matrix for the coefficient functions can be 
expressed in terms of the well-known renormalization constants of QCD and may be found in the papers cited 
above. 

To evaluate the ingredients entering the calculation of T(H — > hadrons) at three-loop level, namely Cf and the 
imaginary part of the correlators {O^O'j), a large amount of diagrams must be calculated. Thus it is necessary 
to use both a generator for the Feynman diagrams and an efficient interface transforming the output to a 
format readable by the integration packages. In the matter at hand one is either confronted with propagator- 
type diagrams or with tadpole integrals where the scales are given by Mh and Mj, respectively. Therefore the 
package GEFIC0M in association with QGRAF (for the generation) and MATAD/MINCER (for the computation) is an 
ideal candidate for the evaluation of the Feynman diagrams contributing to T(H — > hadrons). 





gfffflP 





Figure 21: Typical Feynman diagrams contributing to the coefficients and C%. 
Looped, bold-faced, and dashed lines represent gluons, t quarks, and H bosons, 
respectively. 



In the next section we will give an example of the calculation of three-loop objects in the case of in some 
detail and concentrate in Section 4.4.3| on the evaluation of the massless correlators. GEFICOM was even used 
to compute the 0{ofi s ) contribution to (C^C^) which actually is a four-loop calculation. It is shown how such 



a calculation is nevertheless possible using the three-loop tools provided by GEFICOM. Finally, in Section [4.4.4 
the physical results are briefly presented. 



4.4.2 Coefficient functions 



The top quark dependence of the effective Lagrangian is contained in the coefficient functions Cf. Thus it is 
not surprising that their evaluation can be reduced to massive tadpole integrals with the scale given by M t . 
As an example let us consider C®. There are three independent ways to obtain C®, namely from the sets of 
ggH three-point, gggH four-point, or ggggH five-point diagrams. At three-loop level, these sets contain 657, 
7362, and 95004 diagrams, respectively. A sample diagram of each class is displayed in Fig. |2l] (together with 
a diagram to be evaluated for C®)- The first set requires an expansion up to second order in the external 
momenta, whereas for the second one it suffices to keep only linear terms. For the diagrams of the third class 
all external momenta may be set to zero from the very beginning. However, this advantage is counterbalanced 
by the huge number of contributing diagrams. Note that the depth of the expansion is given by the individual 
terms of the operator 0\ which dictate the structure in the external momenta. 

Ref. [159] deals with the ggH channel, using the package GEFICOM to evaluate the 657 diagrams. The file 
defining the process is rather short because only the package MATAD has to be used: 



*** MATAD 

* power 2 

* gauge 1 

* dalal2 



list = symbolic ; 
lagfile = 'q.lag' ; 
in = h[qlp2] ; 
out = g[ql] , g[q2] ; 
nloop = ; 
options =; 

true = bridge[ g,q,c, 0,0 ]; 



The second line indicates that an expansion in the external momenta up to second order should be performed, 
and dalal2 prompts the application of projectors for the q\-q% terms since we are dealing with on-shell gluons 
(Qi = ?2 = 0)- The fil e containing the propagators and vertices is essentially identical to the pure QCD case 
(cf. Section |3.2.3|) except for the additional coupling of the Higgs boson to top quarks. 

If one chooses a covariant gauge with arbitrary gauge parameter, the CPU time for the calculation is of the 
order of a few days. This is drastically reduced after going, for example, to Feynman gauge. Recently also the 
decay of a pseudo-scalar Higgs boson into gluons was considered [160|. The number of diagrams being the same 
as in the scalar case, there is a slight complication according to the treatment of 7 5 in this case. Therefore the 
7362 diagrams of the three-gluon channel were evaluated in Feynman gauge to have a stringent check of the 



Figure 22: Typical Feynman diagrams contributing to (O^O^)- The solid circles 
represent the operator (re- 



calculation at one's disposal. The CPU time amounts to roughly two weeks. Let us refrain from listing the 



analytical results and refer to [15S, 161] instead. Nevertheless, they are incorporated in the numerical formula 
given in Section 4.4.4 . 



4.4.3 Correlators 



The computation of the one-, two- and three-loop correlators (O^Oj) proceeds in close analogy to the coefficient 
functions described in the previous section. Therefore, let us in this subsection focus on the computation of the 



imaginary part of the four-loop contribution to {0' 2 0' 2 ) which was originally performed in [162]. Some sample 
diagrams are shown in Fig. Although at first sight this seems to require the evaluation of four-loop diagrams, 
the calculation can be reduced to massless propagator-type and massive tadpole integrals at three-loop level. 
The underlying idea is quite simple: since the imaginary part of massless propagator-type diagrams only arises 
from ln(— q 2 ) terms and the logarithms are in one-to-one correspondence to the poles in e it suffices to compute 
the divergent parts of the integrals. It is a feature of dimensional regularization accompanied with MS-like 
schemes that the ultra-violet poles and thus the renormalization constants are polynomials in the masses and 
external momenta. Consequently, for logarithmically divergent diagrams the poles are independent of any mass 
scale. This allows us to assign additional masses to individual lines of the diagrams and to nullify the external 
momentum. If we take one of the lines at the left vertex to be massive and set the external momentum to zero, 
the four-loop integrals can be expanded in terms of three-loop propagator-type diagrams and one-loop tadpole 
integrals for which the technology was described in some detail above. The drawback of the described method is 
that artificial infra-red singularities are introduced which, of course, influence the pole structure. They have to 



be removed using the so-called Infra-red Rearrangement [163] which heavily relies on the R* operation. It was 
developed on a diagram- by-diagram basis in [164]. In order to be able to perform the calculation automatically 
a global version of the R* operation is mandatory [19, 162j , 165]. 



In [162] the computation of Im^^O^) is presented and explicit formulae are given, demonstrating that the 
combination of one-, two- and three- loop massless propagator- type and massive vacuum integrals are adequate 
to get the 0(ag) corrections. The computation was performed with the help of GEFICOM which had to be used 
at four-loop order. A proper flag in the file defining the process tells GEFICOM to introduce the auxiliary mass 
and to rename the momenta in accordance with the factorization taking place. Again we refrain from listing 



explicit expressions and refer to Section 4.4.4 for numerical results. 



4.4.4 Results 



This section summarizes the current knowledge of hadronic Higgs decay by presenting numerical results for the 
decay rates T(H — > gg) and T(H — > bb). 

Note that C\ starts at order a s . Hence the combination C^Im(O^O^) governing the gluonic decay rate of the 
Higgs boson is available up to O(af). Normalized to the Born rate it reads in numerical form (/i = Mjj): 

2 



T(H -> gg) 



1 + 17.917^^1 + 



7T 




156.808 - 5.708 In 



1 + 0.66 + 0.21, 



(64) 



with T Born (H -» gg) = G F M%/{3Giry/2) x (ai 5) (M H )/ir) 2 . In the last line M t = 175 GeV and M H = 100 GeV 
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increases the well-known 0(a s ) enhancement pl],p2| by about one third. If we assume that this trend continues 



to G(cts) and beyond, then Eq. (|64| ) may already be regarded as a useful approximation to the full result. 
Inclusion of the new 0{a 2 s ) correction leads to an increase of the Higgs-boson hadronic width by an amount of 
order 1%. 

Concerning the decay rate into quarks we restrict ourselves to the case of bottom quarks. Inserting numerical 
values into the coefficient functions C\ and C2 and the correlators Iu^C^C^) l!66f | and Im^^O^} 162 , 166| 
leads to: 



A bl \ 1 + 5.667 ag } + 29.147 (ag } ) + 41.758 U% 



(5)V 



(65) 



(ag } ) 2 [3.111 - 0.667 L t ] + (a% ] ) 
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50.474 - 8.167 ^ - 1.278 Lf 



with A b i = 3G F M H ml/ (4:11^/2), L 
and electroweak 
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as ] {M H )/iT. In Eq. ( 
Also mass effects [(351 



|) electromagnetic [(L67 | 
and second order QCD 



158(| corrections have been neglected 
corrections which are suppressed by the top quark mass [34] are not displayed. One observes from Eq. 
that the top- induced corrections at 0(a^) (second line) are of the same order of magnitude as the "massless" 
corrections (first line). 



4.5 Strongly interacting vector bosons 



As an application for the package CompHEP let us consider a scenario where the violation of unitarity in elastic 
scattering of massive vector bosons at high energies is avoided by letting them become strongly interacting. 
This is an alternative approach to the Standard Model solution of introducing a scalar Higgs boson. In [27] 
a new fundamental strong interaction among the vector bosons is assumed that occurs at a scale of order 
1 TeV. The model is described by a Lagrangian with global chiral symmetry which is spontaneously broken 
in order to generate the masses of the vector bosons. The corresponding interaction terms have been worked 
out and implemented in CompHEP. Although this would have been much simpler in unitary gauge due to the 



absence of ghost particles, the authors of [27] decided to use 't Hooft-Feynman gauge as then the size of the 
final expressions is much smaller. One possibility to test the new interaction terms are elastic and quasi-elastic 
2 — > 2 scattering experiments with and Z bosons which requires colliders with a center-of-mass energy in 
the TeV range. Indeed, a future e + e~ linear collider is expected to operate at energies of 1.5 to 2 TeV in a 



second stage [ 169 1 . 

The following processes are considered in 



v e v e ZZ 



e + e~ 

e~e~ — ► v e v e W~W~ 
e + e~ -» e'v e W + Z 
e + e~ — > e + e~ ZZ 



w + w~ -» w + w~ 

ZZ -> ZZ 

W'W~ -» w~w~ 
w+z -» w + z 

ZZ -> ZZ 



(66) 



where after the colon the subprocess characterizing the vector boson scattering is given. However, not only 
those diagrams that contain the desired subprocess contribute. Instead, also background diagrams must be 
computed where, for example, the final state vector bosons are radiated off the fermion lines or one of the 
massive vector bosons is replaced by a photon. Furthermore, diagrams are involved where the neutrinos result 
from the decay of a Z boson which is in turn generated from a pair. In Fig. ^ three sample diagrams are 
listed. 

It is straightforward to go through the menus of CompHEP, to enter the different processes, and to produce 
FORTRAN output for the squared matrix elements. This is then used in the numerical part of CompHEP to 
compute the cross sections. The first run contains no cuts, but computes all matrix elements corresponding to 
the generated set of four- momenta. In subsequent runs, cuts are applied that reduce the background and isolate 
the signal. The data computed in the first run are taken over so that the CPU time is significantly smaller 



in this second step, as already mentioned in Section 3.2.5. For example, the cuts include a lower limit on the 
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Figure 23: Sample diagrams contributing to processes specified in Eq. (|66|) . Dia- 
gram (a) is part of the signal of the vector boson scattering process whereas (b) and 
(c) belong to the background. 
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Figure 24: Distribution of the invariant mass of the produced vector bosons (here 
generically called W) in the process e + e~ — > W + W~u e u e (signal). The cut (shaded 
area) removes events in which the neutrinos are generated through Z decays. 



invariant mass of the v e v e system, or the selection of central events in combination with cuts on the transverse 
momentum of the W^, respectively, the Z boson. Fig. shows how the number of events containing neutrinos 
from Z boson decay (c.f. Fig. p3|(c)) is reduced. Figures of this kind can easily be produced by using the data 
files generated by CompHEP. 

This is not the place to go into details concerning the physical consequences for which we shall refer to [27|. 
The key point, however, is that the couplings of the new interaction terms can be related to the 2^2 vector 
boson scattering amplitudes. Thus, the investigation of the corresponding total cross sections as well as several 
other distributions like the one for invariant masses or transverse momenta provide important tests on the 
mechanism responsible for restoring unitarity at high energies. 



4.6 Higgs mass dependence of two-loop corrections to Ar 

As was already pointed out in the introduction, the dependence of radiative corrections on yet undiscovered 
particles may be exploited in order to gain information on the masses of these particles. This has been 
performed very successfully for the top quark and nowadays, as Mt is measured with reasonable accuracy, the 
same strategy is used to pin down the mass of the Higgs boson. An important quantity in this respect is Ar, 
which comprises the radiative corrections to muon decay. In this way it relates four fundamental parameters of 
the Standard Model to each other — the electromagnetic coupling a, the Fermi constant Gf and the masses 
of the W and Z boson: 

*4-i ) - ^ (1+Ar) 
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Figure 25: Limits on the Higgs mass depending on the values of the top and W 
mass [174|. The value of Ar(M# = 65 GeV) was taken from [173]. 



Gf, a and Mz are very well measured quantities, and the experimental accuracy of My/ is expected to 
be considerably improved within the near future by the current and upcoming measurements at LEP2 and 
TEVATRON. 

Therefore, precise knowledge of the M^-dependence of Ar is highly desirable. The one-loop corrections to 
Ar are known analytically since long |170f| , and two- loop corrections up to next-to- leading order in 1/Mt were 
considered in |171- |173[| . It was the concern of [26 to precisely account for the Mh dependence of the two-loop 



diagrams containing the top quark and the Higgs boson. The idea was to get a result valid for arbitrary values 
of Mi and Mu without performing an expansion in these masses. This obviously leads to two-loop multi-scale 
diagrams, their number being of the order of one hundred. Not only the large number of diagrams but also 
the complicated tensor structure and the evaluation of the scalar two-loop integrals makes the use of highly 
automated software indispensable. A further technical problem arises from the renormalization which has to be 
known up to two loops in the electroweak sector of the Standard Model. In order to investigate the Higgs-mass 
dependece of Ar, in |26]] the subtracted quantity 



Ar(M H ) - Ar(M H = 65 GeV) (68) 
was considered which describes the change in the prediction for Ar when Mh is varied. The value for Ar(M# = 



65 GeV) was taken from [173 



It turned out that the computation of the diagrams contributing to the considered corrections to Ar could be 
reduced to two-loop tadpole integrals and two-point functions. This apparently is a task for the Mathematica 
programs FeynArts and TwoCalc. The various contributions to Ar together with the corresponding counter- 
terms were generated by means of FeynArts. Contraction of Lorentz indices, evaluation of the Dirac algebra 
and the two-loop tensor reduction was performed by TwoCalc. The resulting expressions were split into finite 
and divergent pieces which explicitly demonstrates the cancellation of poles. The scalar one- and two-loop 
integrals were reduced to one-dimensional integrals allowing a fast numerical evaluation to high accuracy with 



the help of special C routines which themselves were fully integrated into the Mathematica environment [102]. 
CPU time for this calculation added up to the order of a few days. 

Instead of displaying the full result, we only note that the resulting prediction for the W mass agrees within a 



few MeV with the one of the heavy-top expansion [173], which is far below the experimental precision at the 



moment. (For a more detailed discussion see [174 and [175].) It is remarkable that the top-induced corrections 



turn out to be maximal at the physical value of the top quark mass within a range of ±50 GeV around this 
physical value. As shown in Fig. ^5], the interrelation between My/ and Mt favours relatively small values 
for Mh-, consistent with other indirect determinations using the radiative shift in the weak mixing angle or 
the leptonic decay rate of the Z boson. Both of these quantities recently have been calculated with similar 
techniques as described in this section [176]. 
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